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Preface 


The field of stochastic computations, in the context of understanding the impact 
of uncertainty on simulation results, is relatively new. However, over the past few 
years, the field has undergone tremendous growth and rapid development. This 
was driven by the pressing need to conduct verification and validation (V&V) and 
uncertainty quantification (UQ) for practical systems and to produce predictions 
for physical systems with high fidelity. More and more researchers with diverse 
backgrounds, ranging from applied engineering to computer science to computa- 
tional mathematics, are stepping into the field because of the relevance of stochas- 
tic computing to their own research. Consequently there is a growing need for an 
entry-level textbook focusing on the fundamental aspects of this kind of stochastic 
computation. And this is precisely what this book does. 

This book is a result of several years of studying stochastic computation and 
the valuable experience of teaching the topic to a group of talented graduate stu- 
dents with diverse backgrounds at Purdue University. The purpose of this book 
is to present in a systematic and coherent way numerical strategies for uncertainty 
quantification and stochastic computing, with a focus on the methods based on gen- 
eralized polynomial chaos (gPC) methodology. The gPC method, an extension of 
the classical polynomial chaos (PC) method developed by Roger Ghanem [45] in 
the 1990s, has become one of the most widely adopted methods, and in many cases 
arguably the only feasible method, for stochastic simulations of complex systems. 
This book intends to examine thoroughly the fundamental aspects of these methods 
and their connections to classical approximation theory and numerical analysis. 

The goal of this book is to collect, in one volume, all the basic ingredients nec- 
essary for the understanding of stochastic methods based on gPC methodology. It 
is intended as an entry-level graduate text, covering the basic concepts from the 
computational mathematics point of view. This book is unique in the fact that it 
is the first book to present, in a thorough and systematic manner, the fundamen- 
tals of gPC-based numerical methods and their connections to classical numerical 
methods, particularly spectral methods. The book is designed as a one-semester 
teaching text. Therefore, the material is self-contained, compact, and focused only 
on the fundamentals. Furthermore, the book does not utilize difficult, complicated 
mathematics, such as measure theory in probability and Sobolev spaces in numer- 
ical analysis. The material is presented with a minimal amount of mathematical 
rigor so that it is accessible to researchers and students in engineering who are 
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xii 

interested in learning and applying the methods. It is the author’s hope that after 
going through this text, readers will feel comfortable with the basics of stochas- 
tic computation and go on to apply the methods to their own problems and pursue 
more advanced topics in this perpetually evolving field. 


West Lafayette, Indiana, USA Dongbin Xiu 

March 2010 
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Chapter One 


Introduction 


The goal of this chapter is to introduce the idea behind stochastic computing in the 
context of uncertainty quantification (UQ). Without using extensive discussions (of 
which there are many), we will use a simple example of a viscous Burgers’ equation 
to illustrate the impact of input uncertainty on the behavior of a physical system and 
the need to incorporate uncertainty from the beginning of the simulation and not as 
an afterthought. 

1.1 STOCHASTIC MODELING AND UNCERTAINTY QUANTIFICATION 

Scientific computing has become the main tool in many fields for understanding the 
physics of complex systems when experimental studies can be lengthy, expensive, 
inflexible, and difficulty to repeat. The ultimate goal of numerical simulations is 
to predict physical events or the behaviors of engineered systems. To this end, ex- 
tensive efforts have been devoted to the development of efficient algorithms whose 
numerical errors are under control and understood. This has been the primary goal 
of numerical analysis, which remains an active research branch. What has been 
considered much less in classical numerical analysis is understanding the impact 
of errors, or uncertainty, in data such as parameter values and initial and boundary 
conditions. 

The goal of UQ is to investigate the impact of such errors in data and subse- 
quently to provide more reliable predictions for practical problems. This topic has 
received an increasing amount of attention in past years, especially in the context 
of complex systems where mathematical models can serve only as simplified and 
reduced representations of the true physics. Although many models have been suc- 
cessful in revealing quantitative connections between predictions and observations, 
their usage is constrained by our ability to assign accurate numerical values to var- 
ious parameters in the governing equations. Uncertainty represents such variability 
in data and is ubiquitous because of our incomplete knowledge of the underlying 
physics and/or inevitable measurement errors. Hence in order to fully understand 
simulation results and subsequently to predict the true physics, it is imperative to 
incorporate uncertainty from the beginning of the simulations and not as an after- 
thought. 

1.1.1 Burgers’ Equation: An Illustrative Example 

Let us consider a viscous Burgers’ equation, 

! u, + uu x = vu xx , x e (—1, 1), 
m(— 1) = 1 , «(!) = -!, 


(1.1) 
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Figure 1.1 Stochastic solutions of Burgers’ equation (1.1) with «(— 1, t) = 1 + 5, where 5 is 
a uniformly distributed random variable in (0, 0.1) and v — 0.05. The solid line 
is the average steady-state solution, with the dotted lines denoting the bounds of 
the random solutions. The dashed line is the standard deviation of the solution. 
(Details are in [123].) 


where u is the solution field and v > 0 is the viscosity. This is a well-known non- 
linear partial differential equation (PDE) for which extensive results exist. The pres- 
ence of viscosity smooths out the shock discontinuity that would develop otherwise. 
Thus, the solution has a transition layer, which is a region of rapid variation and 
extends over a distance of 0(v ) as v 0. The location of the transition layer z, 
defined as the zero of the solution profile u(t,z) = 0, is at zero when the solu- 
tion reaches steady state. If a small amount of (positive) uncertainty exists in the 
value of the left boundary condition (possibly due to some bias measurement or 
estimation errors), i.e., w(— 1) = 1 + <5, where 0 < S <$( 1, then the location of 
the transition can change significantly. For example, if S is a uniformly distributed 
random variable in the range of (0, 0.1), then the average steady-state solution with 
v = 0.05 is the solid line in figure 1.1. It is clear that a small uncertainty of 10 
percent can cause significant changes in the final steady-state solution whose aver- 
age location is approximately at z ^ 0.8, resulting in a 0(1) difference from the 
solution with an idealized boundary condition containing no uncertainty. (Details 
of the computations can be found in [123].) 

The Burgers’ equation example demonstrates that for some problems, especially 
nonlinear ones, a small uncertainty in data may cause nonnegligible changes in the 
system output. Such changes cannot be captured by increasing resolution of the 
classical numerical algorithms if the uncertainty is not incorporated at the begin- 
ning of the computations. 
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1.1.2 Overview of Techniques 

The importance of understanding uncertainty has been realized by many for a long 
time in disciplines such as civil engineering, hydrology, control, etc. Consequently 
many methods have been devised to tackle this issue. Because of the “uncertain” 
nature of the uncertainty, the most dominant approach is to treat data uncertainty as 
random variables or random processes and recast the original deterministic systems 
as stochastic systems. 

We remark that these types of stochastic systems are different from classical 
stochastic differential equations (SDEs) where the random inputs are idealized 
processes such as Wiener processes, Poisson processes, etc., and tools such as sto- 
chastic calculus have been developed extensively and are still under active research. 
(See, for example, [36, 55, 57, 85].) 


1.1. 2.1 Monte Carlo- and Sampling-Based Methods 

One of the most commonly used methods is Monte Carlo sampling (MCS) or one 
of its variants. In MCS, one generates (independent) realizations of random inputs 
based on their prescribed probability distribution. For each realization the data are 
fixed and the problem becomes deterministic. Upon solving the deterministic re- 
alizations of the problem, one collects an ensemble of solutions, i.e., realizations 
of the random solutions. From this ensemble, statistical information can be ex- 
tracted, e.g., mean and variance. Although MCS is straightforward to apply as it 
only requires repetitive executions of deterministic simulations, typically a large 
number of executions are needed, for the solution statistics converge relatively 
slowly. For example, the mean value typically converges as l/\[K , where K is 
the number of realizations (see, for example, [30]). The need for a large number of 
realizations for accurate results can incur an excessive computational burden, espe- 
cially for systems that are already computationally intensive in their deterministic 
settings. 

Techniques have been developed to accelerate convergence of the brute-force 
MCS, e.g., Fatin hypercube sampling (cf. [74, 98]) and quasi Monte Carlo sampling 
(cf. [32, 79, 80]), to name a few. However, additional restrictions are posed based 
on the design of these methods, and their applicability is often limited. 


1.1. 2.2 Perturbation Methods 

The most popular nonsampling methods were perturbation methods, where random 
fields are expanded via Taylor series around their mean and truncated at a certain 
order. Typically, at most second-order expansion is employed because the result- 
ing system of equations becomes extremely complicated beyond the second order. 
This approach has been used extensively in various engineering fields [56, 71, 72]. 
An inherent limitation of perturbation methods is that the magnitude of the uncer- 
tainties, at both the inputs and outputs, cannot be too large (typically less than 10 
percent), and the methods do not perform well otherwise. 
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1.1. 2. 3 Moment Equations 

In this approach one attempts to compute the moments of the random solution 
directly. The unknowns are the moments of the solution, and their equations are 
derived by taking averages of the original stochastic governing equations. For ex- 
ample, the mean field is determined by the mean of the governing equations. The 
difficulty lies in the fact that the derivation of a moment almost always, except on 
some rare occasions, requires information about higher moments. This brings out 
the closure problem, which is often dealt with by utilizing some ad hoc arguments 
about the properties of the higher moments. 

1.1. 2.4 Operator-Based Methods 

These kinds of approaches are based on manipulation of the stochastic operators 
in the governing equations. They include Neumann expansion, which expresses the 
inverse of the stochastic operator in a Neumann series [95, 131], and the weighted 
integral method [23, 24]. Similar to perturbation methods, these operator-based 
methods are also restricted to small uncertainties. Their applicability is often stron- 
gly dependent on the underlying operator and is typically limited to static problems. 

1.1.2. 5 Generalized Polynomial Chaos 

A recently developed method, generalized polynomial chaos (gPC) [120], a gen- 
eralization of classical polynomial chaos [45], has become one of the most widely 
used methods. With gPC, stochastic solutions are expressed as orthogonal poly- 
nomials of the input random parameters, and different types of orthogonal poly- 
nomials can be chosen to achieve better convergence. It is essentially a spectral 
representation in random space and exhibits fast convergence when the solution 
depends smoothly on the random parameters. gPC-based methods will be the focus 
of this book. 

1.1.3 Burgers’ Equation Revisited 

Let us return to the viscous Burgers’ example (1.1), with the same parameter set- 
tings that produced figure 1.1. Let us examine the location of the averaged transition 
layer and the standard deviation of the solution at this location as obtained by differ- 
ent methods. Table 1.1 shows the results by Monte Carlo simulations, and table 1.2 
by a perturbation method at different orders. The converged solutions by gPC (up to 
three significant digits) are obtained by a fourth-order expansion and are tabulated 
for comparison. It can be seen that MCS achieves the same accuracy with O(10 4 ) 
realizations. On the other hand, the computational cost of the fourth-order gPC is 
approximately equivalent to that for five deterministic simulations. The perturba- 
tion methods have a low computational cost similar to that of gPC. However, the 
accuracy of perturbation methods is much less desirable, as shown in table 1.2. In 
fact, by increasing the perturbation orders, no clear convergence can be observed. 
This is caused by the relatively large uncertainty at the output, which can be as high 
as 40 percent, even though the input uncertainty is small. 
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Table 1.1 Mean Location of the Transition Layer (z) and Its Standard Deviation (a z ) by 
Monte Carlo Simulations" 


n = 100 

n = 1,000 

n = 2,000 

n = 5,000 

n = 10,000 

gPC 

z 0.819 

0.814 

0.815 

0.814 

0.814 

0.814 

a z 0.387 

0.418 

0.417 

0.417 

0.414 

0.414 


a n is the number of realizations, S ~ 17(0, 0.1), and v — 0.05. Also shown are the converged 
gPC solutions. 


Table 1.2 Mean Location of the Transition Layer (z) and Its Standard Deviation (cr.) Ob- 
tained by Perturbation Methods" 



k= l 

k = 2 

!*?- 

II 

U> 

Ar = 4 

gPC 

Z 

0.823 

0.824 

0.824 

0.824 

0.814 

O z 

0.349 

0.349 

0.328 

0.328 

0.414 


a k is the order of the perturbation expansion, 8 ~ C7 (0, 0.1), and v = 0.05. Also shown are 
the converged gPC solutions. 

This example demonstrates the accuracy and efficiency of the gPC method. It 
should be remarked that although gPC shows a significant advantage here, the con- 
clusion cannot be trivially generalized to other problems, as the strength and the 
weakness of gPC, or any method for that matter, are problem-dependent. 

1.2 SCOPE AND AUDIENCE 

As a graduate-level text, this book focuses exclusively on the fundamental aspects 
of gPC-based numerical methods, with a detailed exposition of their formulations, 
basic properties, and connections to classical numerical methods. No research top- 
ics are discussed in this book. Although this leaves out many exciting new develop- 
ments in stochastic computing, it helps to keep the book self-contained, compact, 
and more accessible to students who want to learn the basics. The material is also 
chosen and organized in such a way that the book can be finished in a one-semester 
course. Also, the book is not intended to contain a thorough and exhaustive liter- 
ature review. References are limited to those that are more accessible to graduate 
students. 

In chapter 2, we briefly review the basic concepts of probability theory. This is 
followed by a brief review of approximation theory in chapter 3. The material in 
these two chapters is kept at almost an absolute minimum, with only the very ba- 
sic concepts included. The goal of these two chapters is to prepare students for the 
more advanced material in the following chapters. An interesting question is how 
much time the instructor should dedicate to these two chapters. Students taking 
the course usually have some background knowledge of either numerical analy- 
sis (which gives them some preparation in approximation theory) or probability 
theory (or statistics), but rarely do students have both. And a comprehensive cov- 
erage of both topics can easily consume a large portion of class time and leave no 
time for other material. From the author’s personal teaching experience, it is better 
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to go through probability theory rather fast, covering only the basic concepts and 
leaving other concepts as reading assignments. This is reflected in the writing of 
this book, as chapter 2 is quite concise. The approximation theory in chapter 3 de- 
serves more time, as it is closely related to many concepts of gPC in the ensuing 
chapters. 

In chapter 4, the procedure for formulating stochastic systems is presented, and 
an important step, parameterization of random inputs, is discussed in detail. A for- 
mal and systematic exposition of gPC is given in chapter 5, where some of the 
important properties of gPC expansion are presented. Two major numerical ap- 
proaches, stochastic Galerkin and stochastic collocation, are covered in chapters 6 
and 7, respectively. The algorithms are discussed in detail, along with some exam- 
ples for better understanding. Again, only the basics of the algorithms are covered. 
More advanced aspects of the techniques, such as adaptive methods, are left as 
research topics. 

The last chapter, chapter 8, is a slight deviation from the theme of the book be- 
cause the content here is closer to research topics. The topics here, problems in 
random domain, inverse parameter estimation, and “correcting” simulation results 
using data, are important topics and have been studied extensively. The purpose of 
this chapter is to demonstrate the applicability of gPC methods to these problems 
and present unique and effiicient algorithms constructed by using gPC. Neverthe- 
less, this chapter is not required when teaching the course, and readers are advised 
to read it based on their own interests. 


1.3 A SHORT REVIEW OF THE LITERATURE 

Though the focus of this book is on the fundamentals of gPC-based numerical 
methods, it is worthwhile to present a concise review of the notable literature in 
this field. The goal is to give readers a general sense of what the active research 
directions are. Since the field is undergoing rapid development, by no means does 
this section serve as a comprehensive review. Only the notable and earlier work in 
each subfield will be mentioned. Readers, after learning the basics, should devote 
themselves to a more in-depth literature search. 

The term polynomial chaos was coined by Nobert Wiener in 1938 in his work 
studying the decomposition of Gaussian stochastic processes [115]. This was long 
before the phenomenon of chaos in dynamical systems was known. In Wiener’s 
work, Hermite polynomials serve as an orthogonal basis, and the validity of the 
approach was proved in [12]. Beyond the use of Hermite polynomials, the work 
on polynomial chaos referred in this book bears no other resemblance to Wiener’s 
work. In the stochastic computations considered here, the problems we face in- 
volve some practical systems (usually described by partial differential equations) 
with random inputs. The random inputs are usually characterized by a set of ran- 
dom parameters. As a result, many of the elegant mathematical tools in classical 
stochastic analysis, e.g., stochastic calculus, are not directly applicable. And we 
need to design new algorithms that are suitable for such practical systems. 
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The original PC work was started by R. Ghanem and coworkers. Inspired by the 
theory of Wiener-Hermite polynomial chaos, Ghanem employed Hermite polyno- 
mials as an orthogonal basis to represent random processes and applied the tech- 
nique to many practical engineering problems with success (cf. [41, 42, 43, 97]). 
An overview can be found in [45]. 

The use of Hermite polynomials, albeit mathematically sound, presents diffi- 
culties in some applications, particularly in terms of convergence and probability 
approximations for non-Gaussian problems [20, 86]. Consequently, generalized 
polynomial chaos was proposed in [120] to alleviate the difficulty. In gPC, different 
kinds of orthogonal polynomials are chosen as a basis depending on the probability 
distribution of the random inputs. Optimal convergence can be achieved by choos- 
ing the proper basis. In a series of papers, the strength of gPC is demonstrated for 
a variety ofPDEs[119, 121], 

The work on gPC was further generalized by not requiring the basis polynomials 
to be globally smooth. In fact, in principle any set of complete bases can be a 
viable choice. Such generalization includes the piecewise polynomial basis [8, 92], 
the wavelet basis [62, 63], and multielement gPC [110, 11 1]. 

Upon choosing a proper basis, a numerical technique is needed to solve the prob- 
lem. The early work was mostly based on the Galerkin method, which minimizes 
the error of a finite-order gPC expansion by Galerkin projection. This is the stochas- 
tic Galerkin (SG) approach and has been applied since the early work on PC and 
proved to be effective. The Galerkin procedure usually results in a set of coupled 
deterministic equations and requires additional effort to solve. Also, the derivation 
of the resulting equations can be challenging when the governing stochastic equa- 
tions take complicated forms. 

Another numerical approach is the stochastic collocation (SC) method, where 
one repetitively excecutes an established deterministic code on a prescribed node 
in the random space debited by the random inputs. Upon completing the simula- 
tions, one conducts postprocessing to obtain the desired solution properties from 
the solution ensemble. The idea, primarily based on an old technique, the “deter- 
ministic sampling method,” can be found in early works such as [78, 103]. These 
works mostly employed tensor products of one-dimensional nodes (e.g., Gauss 
quadrature). Although tensor product construction makes mathematical analysis 
more accessible (cf. [7]) the total number of nodes grows exponentially fast as 
the number of random parameters grows — the curse of dimensionality. Since each 
node requires a full-scale underlying deterministic simulation, the tensor product 
approach is practical only for low random dimensions, e.g., when the number of 
random parameters is less than 5. 

More recently, there has been a surge of interest in the high-order stochastic 
collocation approach following [118]. A distinct feature of the work in [118] is 
the use of sparse grids from multivariate interpolation analysis. A sparse grid is 
a subset of the full tensor grid and can retain many of the accuracy properties of 
the tensor grid. It can signibcantly reduce the number of nodes in higher random 
dimensions while keeping high-order accuracy. Hence the sparse grid collocation 
method becomes a viable choice in practical simulations. Much more work has 
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since followed, with most focusing on further reduction in the number of nodes 
(cf. [3, 75, 81, 104], 

While most sparse grid collocation methods utilize interpolation theory, another 
practical collocation method is the pseudospectral approach, termed by [116]. This 
approach employs a discrete version of the gPC orthogonal projection operator and 
relies heavily on integration theory. One should keep in mind that in multidimen- 
sional spaces, especially for high dimensions, both interpolation and integration are 
challenging tasks. 

The major challenge in stochastic computations is high dimensionality, i.e.. how 
to deal with a large number of random variables. One approach to alleviate the 
computational cost is to use adaptivity. Current work includes adaptive choice of 
the polynomial basis [33, 107], adaptive element selection in multielement gPC 
[31, 110], and adaptive sparse grid collocation [35, 75, 104]. 

Applications of these numerical methods take a wide range, a manifestation of 
the relevance of stochastic simulation and uncertainty quantification. Here 
we mention some of the more representative (and published) work. It includes 
Burgers’ equation [53, 123], fluid dynamics [58, 61, 64, 70, 121], flow-structure 
interactions [125], hyperbolic problems [17, 48, 68], material deformation [1, 2], 
natural convection [35], Bayesian analysis for inverse problems [76, 77, 112], 
multibody dynamics [89, 90], biological problems [37, 128], acoustic and electro- 
magnetic scattering [15, 16, 126], multiscale computations [5, 94, 117, 124, 129], 
model construction and reduction [26, 40, 44], random domains with rough bound- 
aries [14, 69, 102, 130], etc. 


Chapter Two 


Basic Concepts of Probability Theory 


In this chapter we collect some basic facts needed for stochastic computations. 
Most parts of this chapter can be skipped, if one has some basic knowledge of 
probability theory and stochastic processes. 


2.1 RANDOM VARIABLES 

The outcome of an experiment, a game, or an event is random. A simple example 
is coin tossing: the possible outcomes, heads or tails, are not predictable in the 
sense that they appear according to a random mechanism that is too complex to be 
understood. A more complicated experiment is the stock market. There the random 
outcomes of brokers’ activities, which in fact represent the economic environment, 
political interests, market sentiment, etc., are reflected by share prices and exchange 
rates. 

The mathematical treatment of such kinds of random experiments requires that 
we assign a number to each random outcome. For example, when tossing a coin, 
we can assign 1 for heads and 0 for tails. Thus, we obtain a random variable 
X = X(w) e {0, 1}, where a> belongs to the outcome space £T2 = {heads, tails}. 
In the example of the stock market, the value of a share price of a stock is already a 
random variable. The numbers X(a>) provide us with information about the exper- 
iment even if we do not know precisely what mechanism drives the experiment. 

More precisely, let !T2 be an abstract space containing all possible outcomes o> 
of the underlying experiment; then the random variable X — X(co) is a real- valued 
function defined on JT2 . Note here that for the abstract space Q, it does not really 
matter what the ca are. 

To study problems associated with the random variable X, one first collects rele- 
vant subsets of Q, the events, in a class T called a a -field or a -algebra. In order for 
T to contain all relevant events, it is natural to include all the co in the event space 

and also the union, difference, and intersection of any events in T , the set Q, and 
its complement, the empty set 0. 

If we consider a share price A of a stock, not only the events {u> : X(ai) = c} 
should belong to T but also 

{co : a < X(a>) < b], {co : b < X{co)}, {co : X(co) < a], 

and many more events that can be relevant. And it is natural to require that elemen- 
tary operations such as fT (J, c on the events of T will not land outside the class 
T . This is the intuitive meaning of a o -field T . 
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Definition 2.1. A a -field T (on £2) is a collection of subsets of £2 satisfying the 
following conditions: 

• It is not empty: 0 e T and £2 e 

• If A e T, then A c e T . 

• If A\, A 2 , ■ ■ ■ , e T, then 

OO OO 

Aj e T and (^| A t e T . 

1=1 /= 1 

Example 2.2 (Some elementary a -fields). The following collections of subsets of 
£2 are cr-fields: 


T x = {0, £2} 

•F 2 = {0, Q,A,A c }, where ^ ^ 0, A ^ £2, 

T^=2 q = [A : A C £2}. 

T\ is the smallest cr-field on £2, and is the biggest one containing all possible 
subsets of £2 and is called the power set of £2. 

In practice, the power set is in general too big. One can prove that, for a given 
collection C of subsets of £2, there exists a smallest cr-field cr(C) on £2 containing 
C. We call cr(C) the a -field generated by C. 

Example 2.3. In example 2.2, one can prove that 

Tx =<T({0}), T 2 = a({A}), J~3 = (j(J-f). 


2.2 PROBABILITY AND DISTRIBUTION 

The concept of probability is used to measure the likelihood of the occurrence 
of certain events. For example, for the fair coin toss described in the previous 
section, we assign the probability 0.5 to both events, heads and tails. That is, 
P({a> : X(u>) = 0}) = P({a> : X(u>) = 1}) = 0.5. This assignment is based 
on empirical evidence: if we flip a fair coin a large number of times, we expect that 
about 50 percent of the outcomes will be heads and about 50 percent will be tails. 
In probability theory, the law of large numbers gives the theoretical justification for 
this empirical observation. 

Some elementary properties of probability measures are easily summarized. For 
events A, B e T, 

P(A U B)= P(A) + P(B ) - P(A n B), 
and if A and B are disjoint, 

P(AUB)= P(A) + P(B). 

Moreover, 


P(A C ) = 1 - P(A), P(£ 2) = 1, P(0) = 0. 
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Definition 2.4 (Probability space). A probability space is a triplet (£2, IF, P ) 
where £2 is a countable event space, T C is the o -field ofQ, and P is a proba- 
bility measure such that 

1. 0 < P(A) < 1 ,VA e T. 

2. P(Q) = 1. 

3. For A\, A 2 , ■ ■ . € T and A ,■ fl Aj = 0, V/ ^ j, 

( OO \ OO 

U 4J = '£p( a i ). 

Definition 2.5 (Distribution function). The collection of the probabilities 

F x (x) = P(X <x)= P({a> : X(ea) < x}), ieI, (2.1) 

is the distribution function F x of X. 

It yields the probability that X belongs to an interval (a, b\. That is, 

P({a > : a < X(co) < b}) = F x (b) - F x (a), a <b. 

Moreover, we obtain the probability that X is equal to a number 
P(X = x)= F x (x) - lim F x ( x - e). 

e—>0 

With these probabilities we can approximate the probability of the event {o> : 
X(a>) e B } for very complicated subsets B of K. 

Definition 2.6 (Distribution). The collection of the probabilities 
P X (B ) = P(X e B) = P({w : X(co e B}) 
for suitable subsets B C® is the distribution of X. 

The suitable subsets of R are called Bore l sets. They are sets from B — cr({(a, b] : 
— 00 < a < b < 00 }), the Borel a -field. 

The distribution P X and the distribution function F x are equivalent notions in 
the sense that both of them can be used to calculate the probability of any event 
{XeB}. 

2.2.1 Discrete Distribution 

A distribution function can have jumps. That is, 

F x {x) = ^2 Pk, x e R, (2.2) 

k:x/c<x 

where 

OO 

0 < Pk < 1 , VA:, y jPk =l. 

k= 1 

The distribution function (2.2) and the corresponding distribution are discrete. A 
random variable with such a distribution function is a discrete random variable. 

A discrete random variable assumes only a finite or countably infinite number of 
values xi,X 2 , . ■ ■ and with probability — P(X — xfi). 
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Figure 2.1 Probability distribution functions. Left: binomial distribution with n = 10, p = 
0.5. Right: Poisson distribution with X = 3. 


Example 2.7 (Two important discrete distributions). Important discrete distri- 
butions include the binomial distribution Bin, p) with parameters n e No = 
{0, 1, . . . } and p e (0, 1): 

P(X = k)= ff)/(l -/>)"-*, At = 0, 1, ...,«, 

and the Poisson distribution P(F) with parameter '/, > 0: 

P(X= k) = e~ 1 — , it =0,1 

it! 

Graphical illustrations of these two probability distributions are shown in figure 2. 1 . 


2.2.2 Continuous Distribution 

In contrast to discrete distributions and random variables, the distribution function 
of a continuous random variable does not have jumps; hence 

P(X=x) = 0, VreR, 


or, equivalently. 


lim F x (x + e) = F x (x), Vx; (2.3) 

a — ^0 

i.e., a continuous random variable assumes any particular value with probability 0. 
Most continuous distributions have a density f x : 

Fx(x) = f fx(y)dy , ret, (2.4) 

J — oo 


where 


fx(x) >0, Vx e K, 



fx(y)dy = 1 . 


Example 2.8 (Normal and uniform distributions). An important continuous dis- 
tribution is the normal or Gaussian distribution Af (ft, cr 2 ) with parameters /iet, 
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Figure 2.2 Normal distribution function with p = 0, a 2 = 1. 


a 2 >0. Its density is 

fx(x) = 


/in a 2 


exp 


( X - p) 1 

2(7 2 


The density of A/"(0, 1) is shown in figure 2.2. 

The uniform distribution U(a, b) on (a, b ) has density 


fx(x) 


which is a constant inside (a, b). 


x e (a, b), 
0, otherwise, 


(2.5) 


2.2.3 Expectations and Moments 


Important characteristics of a random variable X include its expectation, variance, 
and moments. The expectation or mean value of a random variable X with density 
fx is 


px — E[X| 



The variance of X is defined as 

o\ = var(X) = 

J —OO 

The m\h moment of X for m e N is 


/ OO 

(x - p x ) 2 fx{x)dx. 

-C 

S 

/ OO 

x m f x (x)dx. 

-OO 


For a real- valued function g, the expectation of g( X) is 

poo 

E[g(JQ] = / g(x)f x (x)dx. 
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Similarly, for a discrete random variable X with probabilities pk = P(X = xf), 
we have 

OO 

fix = MX\ = ^2x k p k , 
k = 1 
oo 

o\ = var(Z) = ^2(x k - iix) 2 Pk , 

*=i 

oo 

£=1 

oo 

E[*(*)] = £>(**)/%■ 

*=i 

Often we study E [(X — Hx)"' )- the centered moments. The expectation or mean 
/r x, is often regarded as the “center” of the random variable X or the most likely 
value of X. The variance a x , or more precisely the standard deviation n x , de- 
scribes the spread or dispersion of the random variable X around its mean fix- It is 
straightforward to show that 

ct 2 =E[(X- fx x ) 2 ] 

= E[X 2 - 2/j.xX + fi 2 x \ = E[X 2 ] - 2pr x + ii 2 x 
= E[X 2 ]-,i 2 x . 

Often this is memorized as “variance equals the mean of the square minus the 
square of the mean.” 

Example 2.9 (Moments of a Gaussian random variable). If a random variable X 
has a normal distribution with density (2.5), then we have 

lix = E[Z] — p,, 

a x = var(JT) = cr, 

E [(X — /x) 2 " -1 ] = 0, n = 1,2,..., 

E \{X — /x) 2 "] = 1 - 3 ■ 5 - - ■ (2« — 1) -a 2n , n = 1,2, .... 

Therefore, the mean and variance of a normal distribution can completely charac- 
terize all of its moments. 

2.2.4 Moment-Generating Function 

Definition 2.10. The moment-generating function for a random variable X(o>) is 
defined as nix(t) = E[e' A ], It exists if there exists b > 0 such that mx{t) is finite 
for |f | < b. The reason mx(t ) is called a moment-generating function is because 

d k m x {t) 


b-k = 


£ = 0 , 1 ,..., 
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where pt k — E[T*] is the kth moment of X. The relationship can be seen as follows: 

mx(t) = E[e ,x ~\ = J e tx px(x)dx 

f °° tx k 00 1 r 

= 22 ~rp Px(x)dx = — / t k x k p x (x)dx 

k= 0 ' k=0 K ' J 

“ t k p k t 2 

= y, —rr~ = do + tp.\ + — fi 2 -\ ■ (2.6) 

z — ' k\ 2 

k = 0 

If X ~ A/"(0, cr 2 ) w a Gaussian random variable, then mx(t ) = e !T " , ~/ 2 . 

2.2.5 Random Number Generation 

One of the basis tasks in stochastic simulations is to generate a sequence of ran- 
dom numbers satisfying a desired probability distribution. To this end, one first 
seeks to generate a random sequence with a common uniform distribution in (0, 1). 
There are many available algorithms that have been well studied. The algorithms in 
practical implementations are all deterministic (typically using recursion) and can 
therefore at best mimic properties of uniform random variables. For this reason, 
the sequence of outputs is called a sequence of pseudorandom numbers. Despite 
some defects in the early work, the algorithms for generating pseudorandom num- 
bers have been much improved. The readers of this book will do well with existing 
software, which is fast and certainly faster than self-made high-level-language rou- 
tines, and will seldom be able to improve it substantially. Therefore, we will not 
spend time on this subject, and we refer interested readers to references such as 
[38, 59, 65, 87], 

For nonuniform random variables, the most straightforward technique is via the 
inversion of a distribution function. Let Fx(x) = P(X < x) be the distribution 
function of X. For the simple case where F x is strictly increasing and continuous, 
then x = Ff l ( u ) is the unique solution of F x (x) — u, 0 < u < 1. For distributions 
with nonconnected support or jumps, F x is not strictly increasing, and more care is 
needed to find its inverse. We choose the left-continuous version 

Ff\u) = inf}* : F x (x) > u}. (2.7) 

We state here the following results that justify the inversion method for generating 
nonuniform random numbers. 

Proposition 2.11. Let F x (x) = P(X < x) be the distribution function of X. Then 
the following results hold. 

• u < F x (x) Ff l (u) < x. 

• IfU is uniform in (0, 1), then Ff l (U) has distribution function F x . 

• If Fx is continuous, then F x (x) is uniform in (0, 1). 

The part that is mainly used in simulation is the second statement, which allows 
us to generate X as Fx 1 (U). The most common case is an F x that is continuous 
and strictly increasing on an interval. 
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Example 2.12 (Exponential random variable). Let Abe a random variable with 
exponential distribution whose probability density is f x {x ) = ae~ ax . Here a > 0 
is a parameter that is often called the rate. It is easy to see that E[A] = 1 /a, the 
inverse of the rate. Then the inverse of F x is F x ( u ) = — log(l — x)/a. So we 
can generate X by inversion: X = — log(l — U)/a. In practice, one often uses the 
equivalent X = — log(U)/a. 

A main limitation of inversion is that quite often Fx is not available in explicit 
form, for example, when A is a normal random variable. Sometimes approxima- 
tions are used. 


Example 2.13 (Approximate inversion of normal distribution). Let A be a 

Gaussian random variable with zero mean and unit variance; i.e., X ~ AGO. 1). 

Its probability density function is f x (x) = -k=e~ x 7 2 , and there is no explicit for- 

V2 n 

mula for F x (x). The following two approximations to F x ! are quite simple and 
accurate. 


F x 1 (u) & sign(t« — 1 /2) 


Co + C\t + C 2 f 2 \ 
1 T d\t -f- dit~ + d\ t ^ J 


(2.8) 


where t = (— ln[min(M, (1 — w))] 2 ) 1 / 2 and co = 2.515517, ci = 0.802853, d\ = 
1.432788,(72 = 0.189269,(73 = 0.001308. The formula has absolute error less than 
4.5 x 10~ 4 ([50]). Or, 


F x (u) ^ y + 


po + piy + pit + piy 


P4 y 


qo + qiy + qiy- + 93 / + q^r 


0.5 < u < 1, 


(2.9) 


where v = v /— 2 log( I — u ) , the case of 0 < u < 0.5 is handled by symmetry, and 
Pk , qk are given in the accompanying table. 


k 

Pk 

qk 

0 

-0.322232431088 

0.099348462606 

1 

-1 

0.588581570495 

2 

-0.342242088547 

0.531103462366 

3 

-0.0204231210245 

0.10353775285 

4 

-0.0000453642210148 

0.0038560700634 


Other techniques exist for the generation of nonuniform random variables, most 
notably acceptance-rejection algorithms. We will not engage in more in-depth dis- 
cussions here and refer interested readers to comprehensive references such as 
[25, 38, 52], 


2.3 RANDOM VECTORS 

We say X = (Aj , . .., X n ) is an n -dimensional random vector if its compo- 
nents X\ , , X n are one-dimensional real-valued random variables. Therefore, a 
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random vector is nothing but a collection of a finite number of random variables. 
Similarly, we can also define concepts such as distribution function, moments, etc. 

Definition 2.14. The collection of the probabilities 

F x (x) = P(Xi <x u ...,X n <x„), x=(xi eR", (2.10) 


is the distribution function Fx ofX. 

If the distribution of a random vector X has a density fx, we can represent the 
distribution function F\ as 


Fx(*i, . . . ,x„) = 


-r-r 


fxiyu 


. }’n)dy i ■ • • dv„, 


where the density is a function satisfying 

/x(x) >0, Vxe 

and 


/ OO PC 

-OO J —C 


fx(y \ , ■ ■ • , y„)dyi ■■■dy„ = 1 . 


If a vector X has density fx, then all of its components X,. the vectors of the 
pairs ( X , , Xf), triples (X , , Xj, Xf), etc., have a density. They are called marginal 
densities. For example, 


fxfxi) = 


/ OO PC 

■■■/ 

-oo J — c 


fxiyu ■ ■ ■, y n )dy\ ■ ■ ■ dyi_idy i+ \ ■ ■ • dy n 


The important statistical quantities of a random vector include its expectation, 
variance, and covariance. The expectation or mean value of a random vector X is 
given by 

Ai X = E[X] = (E[X!],...,E[X„]). 

The covariance matrix of X is defined as 


C x = (co v(X t ,Xj))l J=1 , (2.11) 

where 


co v(X M Xj) = E[(*i- - Hx,)iXj - Hxf] = K{XiXj) - p Xi p Xj , (2.12) 

is the covariance of X, and X r Note that cov( X,, X, ) = a\.. 

It is also convenient to standardize covariances by dividing the random variables 
by their standard deviations. The resulting quantity 

cov(2fi, Xf) 

corr(Xi, XX) = v (2.13) 

V Xl VX 2 

is the correlation coefficient. An immediate fact following the Cauchy-Schwarz 
inequality is 

-1 < corr(Au, Xf) < 1. (2.14) 

We say the two random variables are uncorrelated if corr(X! , X 2 ) = 0, and strongly 
correlated if |corr(Ai, Xf)\ & 1. 
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Example 2.15 (Gaussian random vector). A Gaussian or normal random vector 
has a Gaussian or normal distribution. The n -dimensional Gaussian distribution is 
defined by its density 



M ' )= (2„ )-/>(!, Cx)''» exp 


where /i\ e E" is the expectation of X and Cx is the covariance matrix. Thus, the 
density of a Gaussian vector (hence its distribution) is completely determined via its 
expectation and covariance matrix. If Cx = I„, the n -dimensional identity matrix, 
the components are called uncorrelated and the density becomes the product of n 
normal densities: 


/x(x) = f Xl (xi) ■ ■ ■ fx n (x n ) , 


where fx t (Xi ) is the normal density of Af( /x x, ■ cr| ). An important and appealing 
property of Gaussian random vectors is that they remain Gaussian under linear 
transformation. 

Theorem 2.16. Let X = (X \, . . . , X n ) be a Gaussian random vector with distribu- 
tion N(p, C) and let A be an m x n matrix. Then AX r has an AT(A/j. t , ACA r ) 
distribution. 

The proof is left as an exercise. 


2.4 DEPENDENCE AND CONDITIONAL EXPECTATION 

Intuitively, two random events are called independent if the outcome of one event 
does not influence the outcome of the other. More precisely, we state the following. 

Definition 2.17. Two events A\ and /1 2 are independent if 


P{A l C\A 2 ) = P{A x )P(A 2 ). 


Definition 2.18. Two random variables X \ and X 2 are independent if 

P(Xi eB u X 2 e B 2 ) = P(X 1 e B { )P(X 2 e B 2 ) 

for all suitable subsets B\ and B 2 o/R. This means that the events {X\ e B\) and 
{ X 2 e ZR | are independent. 

Alternatively, one can define independence via distribution functions and densities. 
The random variables X\, . . . , X n are independent if and only if their joint distri- 
bution function can be written as 

F Xl ,....x n {xi, ■ ■ -,x„) = F Xl (xi ) ■ ■■ F Xn (x„), (xi,...,x n ) e R". 

If the random vector X = (X\ , . . . , X n ) has density /x, then X\, . . . , X n are inde- 
pendent if and only if 


f Xl ,...,x„ (xi,...,x n ) = fxfxi) ■ ■ ■ fx„{x n ), (xi,...,x„) eR". (2.16) 


CONCEPTS OF PROBABILITY THEORY 


19 


It follows immediately that if X \ , . . . , X n are independent, then for any real- valued 
functions gi, . . . , g n , 

E[gi(Xi) ■ ■ ■ g n (X n )] = E[a,(X,)] ■ ■■E[g„(X„)], 

provided the considered expectations are all well defined. Hence, if X\ and X 2 are 
independent, then 

corr(Xi, X 2 ) = cov(A), X 2 ) = 0. 

This implies that independent random variables are uncorrelated. The converse is 
in general not true. 

Example 2.19. Let X be a standard normal random variable. Since X is symmetric 
(its density is an even function), so is X 3 . Therefore, both X and X 3 have expecta- 
tion zero. Thus, 

cov(X, X 2 ) = E[X 3 ] - E[JT]E[X 2 ] = 0, 

which implies that X and X 2 are uncorrelated. However, they are clearly dependent, 
in the sense that knowing X determines X 2 completely. For example, since {X e 
[-1, 1]} = {X 2 e [0, 1]}, we have 

P(Xe [-1, 1],X 2 e [0, l]) = P(Xe [-1, 1]) 

> P(Xe [-1, 1 })P(X 2 e [0, 1]) 

= (P(Xe [-1, l])) 2 . 

Example 2.20. Let X be an n -dimensional Gaussian random vector with density 
(2.15). The components are uncorrelated when corr( X, . Xj) = cov(A), Xf) = 0 
for i ^ j. This means that the correlation matrix is diagonal. In this case, the den- 
sity of X can be written in the product form of (2. 16), and therefore the components 
are independent. Thus, uncorrelation and independence are equivalent notions for 
Gaussian distributions. 


Another concept describing relations among random events or random variables 
is the conditional expectation. This is an extremely important subject in probability 
theory, though it is not of as much significance in this book. 

From elementary probability theory, the conditional probability of A given B is 


P(A\B) = 


P(A rr B) 

~p(b) 


Clearly, 


P( A | B) = P( A ) if and only if A and B are independent. 


Given that P( B) > 0, we can define the conditional distribution function of a 
random variable X given B, 

P(X < x, B) 


Fx(x\B) = 


P(B) 


and also the conditional expectation of X given B, 

E [XI B ] 


E[X\ B] = 


P{B) 


(2.17) 
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where 


Ib((o) = 


co e B, 
co £ B, 


denotes the indicator function of the event B. For the moment let us assume that 
£2 = R. If X is a discrete random variable, then (2.17) becomes 


^ P({« : X(oo) = x k ] n B) ^ Im 

E[X\B] = 2_^x k — — = 2_x k P(X = x k \B). 


P(B) 

If X has density f x , then (2.17) becomes 
1 


E[X\ B] = 


P(B ) 


/ 


xI B (x)fx(x)dx = 


-M 


xfx(x)dx. 


2.5 STOCHASTIC PROCESSES 

In many physical systems, randomness varies either continuously or discretely over 
physical space and/or time. Therefore, it is necessary to study the distribution and 
evolution of the random variables that describe the randomness as functions of 
space and/or time. A mathematical model for describing this is called a stochastic 
process or a random process. 

A stochastic process is a collection of random variables 

(X t , t eT) = (Xfico), t e T, co e ft) 

defined on some space Here t is the index of the random variable X. The index 
set T can be an interval, e.g., T = [a, b], [a, b) or [a, oo) for a < b. Then X is 
a continuous process. If T is a finite or countably infinite set, then A is a discrete 
process. Very often the index t of X t is referred to as time. However, one should 
keep in mind that it is merely an index and can be a space location as well. 

A stochastic process X can be considered as a function of two variables. 

• For a fixed index t, it is a random variable: 

Xf = Xf ((d) , co G £2 . 

• For a fixed random outcome oo e £2, it is a function of the index (time): 

X, = Xfoo), teT. 

This is called a realization, a trajectory, or a sample path of the process X t . 

It is then natural to seek the statistics of the stochastic process. The task is more 
complicated than that for a random vector. For example, a process X, with an in- 
finite index set T is an infinite-dimensional object. Mathematical care is needed 
for such objects. A simpler approach, which suits practical needs well, is to in- 
terpret the process as a collection of random vectors. In this way, we study the 
finite-dimensional distributions of the stochastic process X t , which are defined as 
the distributions of the finite-dimensional vectors 

{x h ,...,x tn ). 


t\ , . . . , t n g r, 
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for all possible choices of the index t\, ... ,t n e T and every n > 1 . These are 
easier to study and indeed determine the distribution of X,. In this sense, we refer 
to a collection of finite-dimensional distributions as the distribution of the stochastic 
process. 

A stochastic process X = (X t , t e T) can be considered a collection of random 
vectors (X h , .X, n ) for t\, ... ,t n e T and n > 1 . We can then extend the defini- 
tions for expectation and covariance matrices for the random vector to the process 
and consider these quantities as functions of t e T . The expectation function of X 
is given by 


p-x(t) = p.x, — EtA’t], t e T. 


The covariance function of X is given by 

C x (t, s ) = cov(X t , X s ) = E[(A, - pt x {t)){X s - ti x (s))l t, s e T. 

The variance of X is given by 

o 2 x (t) = C x (t , t) = var(X,), t e T. 

These are obviously deterministic quantities. The expectation function p. x {t) is a 
deterministic path around which the sample paths of X are concentrated. Note that 
in many cases p, x (t) may not be a realizable sample path. The variance function 
can be considered a measure of the spread of the sample paths around the expecta- 
tion. The covariance function is a measure of dependence on the process X. 

The process X = ( X t , t e T) is called strictly stationary if the finite-dimensional 
distributions of the process are invariant under shifts of the index t: 

(X h , . . . , X t J = (X h+h , . . . , X tn+h ) 

for all possible choices of indices t\, ... ,t n e T, n > 1 and h such that t\ + 
h, ... , f„ + h e T . Here = stands for the identity of the distributions; see the 
following section for the definition. In practice, a weaker version of stationarity is 
often adopted. A process X is called stationary in the wide sense or second-order 
stationary if its expectation is a constant and covariance function C x(t . s) depends 
only on the distance 1 1 — ,s T | . For a Gaussian process, since its mean and covariance 
function can fully characterize the distribution of the process, the two concepts for 
stationarity become equivalent. 

A large class of (extremely) useful processes can be constructed by imposing a 
stationary (strictly or in the wide sense) condition on the increment of the processes. 
Examples include the homogeneous Poisson process and Brownian motion. Exten- 
sive mathematical analysis has been conducted on such processes by using nonele- 
mentary facts from measure theory and functional analysis. This, however, is not 
the focus of this book. We refer interested readers to the many excellent books such 
as [36, 55], 
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2.6 MODES OF CONVERGENCE 

We now introduce concepts of main modes of convergence for a sequence of ran- 
dom variables X\ , Xi , . . . . 

Definition 2.21 (Convergence in distribution). The sequence { X n ) converges in 
distribution or converges weakly to the random variable X, written as X n —>■ X, if 
for all bounded and continuous functions f, 

E [/(*„)] -* E(/(W)], n -* 00. 

Note that X n —>■ X holds if and only if for all continuous points x of the distribution 
function Fx the relation 


F Xn (x ) F x (x), n 00 , 

is satisfied. If Fx is continuous, this can be strengthened to uniform convergence: 
sup \ F Xn (x) — F x (x)\ — > 0 , n — >■ 00. 

We state the following useful theorem without proving it. 

Theorem 2.22. Let X n and X be random variables with moment-generating func- 
tions mx„{t) and mxif), respectively. If 

lim mx.(t ) = mx(t ), Vt, 

w— s>-oo 

then X n —*■ X as n — y 00. 

Definition 2.23 (Convergence in probability). The sequence { X „ ) converges in 

p 

probability to X, written as X n — > X, if for all positive e, 

P( \X„ - X\) > e) -> 0, n -» 00. 

Convergence in probability implies convergence in distribution. The converse is 
true if and only if X — x for some constant x. 

Definition 2.24 (Almost sure convergence). The sequence \ X„} converges almost 
surely (a.s.), or with probability 7, to the random variable X, written as X n X, 
if the set of co with 

X n (co) -» X(a>), n -» 00, 


has probability 1. 

This implies that 

P(X n -> X) = P({co : X n (oj) -> X(m)}) = 1 . 

Convergence with probability 1 implies convergence in probability, hence conver- 
gence in distribution. Convergence in probability does not imply a.s. convergence. 

P as 

However, X n -> X implies that X„ t — » A for a suitable subsequence { X„ t j . 
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Definition 2.25 ( L p convergence). Let p > 0. The sequence { X n ) converges in 

L p 

L p , or in the pth mean, to X, written as X n -* X, ifE[\X n \ p + |X| P ] < oo for all 
n and 

E[|A„ - X\ p ] -» 0, n — y oo. 

The well-known Markov’s inequality ensures that P( \X„ — X\ > e) < e _ CE[| X n — 

L p P 

X\ p ] for positive p and e. Thus, X n — >• X implies X n — > X. The converse is in 
general not true. 

For p = 2, we say that X n converges to X in mean square. Mean-square conver- 
gence is convergence in the Hilbert space 

L 2 = L 2 (Q, T, P) = {X: E[X 2 ] < oo} 

endowed with the inner product (X, Y) = E[XT] and the norm || X\\ = *J(X, X). 

Convergence in distribution and convergence in probability are often referred 
to as weak convergence, whereas a.s. convergence and L p convergence are often 
referred to as strong convergence. 


2.7 CENTRAL LIMIT THEOREM 


The celebrated central limit theorem (CLT) plays a central role in many aspects of 
probability analysis. Here we state a version of the CLT that suits the exposition of 
this book. 


Theorem 2.26. Let X\, X 2 , ...,X„ be independent and identically distributed 
(i.i.d.) random variables with E[A}] = p and var(Xj) — o 2 < 00 . Let 


x=-Yx t 

1=1 


and let 


U n = J7t 


X- 


Then the distribution function of U n converges to an A/”(0, 1) distribution function 
as n 00 . 


Proof For all i, let Z, = Xi a ^ ; then E[Z,] = 0 and var(Z, ) = 1. Also, let = 
E[Z ( 3 ] for all i. Then, by following the property of moment-generating function 
(2.6), we have 

t 2 / 3 

=l + - + —p 3 H . 

By definition, 


U„ = Vn 


X- p 


1 Xj — np \ _ z 

yfn \ o ) yfn “ 
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We immediately have 



Letting z = £ + 3 ^ 372/^3 H and using ln(l + z) = z- § + ^ 

have 


In {ma n ( 0 ) = n 




t 3 H 3 

3 !« 3 / 2 


It is obvious that 



lim m Un (t) - — , 

w — >00 2 

which is the moment-generating function of a unit Gaussian random variable 
A/"(0, 1). The theorem is then proved by virtue of theorem 2.22. 


This immediately implies that the numerical average of a set of i.i.d. random vari- 
ables {X ,}" =1 will converge, as n is increased, to a Gaussian distribution 
Af(n,<j 2 /n), where /i and a 2 are the mean and variance of the i.i.d random 
variables. 


Chapter Three 


Survey of Orthogonal Polynomials and 
Approximation Theory 


In this chapter we review the basic aspects of orthogonal polynomials and approxi- 
mation theory. We will focus exclusively on the univariate case, that is, polynomials 
and approximations on the real line, to establish the fundamentals of the theories. 


3.1 ORTHOGONAL POLYNOMIALS 

We first review the basics of orthogonal polynomials, which play a central role in 
modern approximation theory. The material is kept to a minimum to satisfy the 
needs of this book. More in-depth discussions of the properties of orthogonal poly- 
nomials can be found in many standard books such as [10, 19, 100]. 

From here on we adopt the standard notation by letting N be the set of positive 
integers and letting No be the set of nonnegative integers. We also let AT = No = 
{0, 1 , . . . } or AT = {0, 1 , . . . , N} be an index set for a finite nonnegative integer N. 


3.1.1 Orthogonality Relations 

A general polynomial of degree n takes the form 

Q n (x) = a„x n + + • • • + a\x + ciq, a n ^ 0, (3.1) 

where a n is the leading coefficient of the polynomial. We denote by 

P, W = + +... + £!, + fi 

a„ a„ a n a n 

the monic version of this polynomial, i.e., the one with the leading coefficients 
equal to 1 . 

A system of polynomials {Q„(x),n e A f} is an orthogonal system of polyno- 
mials with respect to some real positive measure a if the following orthogonality 
relations hold: 


L 


Q n (x)Q m (x)da(x) = y n &„ 


,neAf, 


(3.2) 


where S m „ = 0 if m ^ n and S mn = 1 if m — n, is the Kronecker delta func- 
tion, S is the support of the measure a, and y„ are positive constants often termed 
normalization constants. Obviously, 


l 


Yn = / Q,,(x)da(x), 


e Af. 
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If y n = 1, the system is orthonormal. Note that by defining Q n (x ) = Q„(x)/ ' JYn, 
the system { Q n } is orthonormal. 

The measure a usually has a density w(x) or is a discrete measure with weight 
Wi at the points x,. The relations (3.2) then become 


L 


Qn (x) Q m (x)u)(x)dx = Yn&n 
in the former case and 

Y Qn(Xi)Q m (Xi)Wi = YnSmn 


, n e J\f, 




in the latter case where it is possible that the summation is an infinite one. 
If we define a weighted inner product 


(u, v)da = I u(x)v(x)da(x), 

Js 

which in continuous cases takes the form 

(u, v) w = / u(x)v(x)w(x)dx 

Js 

and in discrete cases takes the form 

(u, v) w = Y, u(xj)v(xi)Wj, 

then the orthogonality relations can be written as 

(Q m Q n ) w = Yn&mn m,H € A/", 

where 

Yn={Qn,Qn) W =\\Qn\\l, « £ AT. 


(3.3) 

(3.4) 

(3.5) 

(3.6) 

(3.7) 

(3.8) 

(3.9) 


3.1.2 Three-Term Recurrence Relation 

It is well known that all orthogonal polynomials { Q„(x)} on the real line satisfy a 
three-term recurrence relation 

-xQ n (x) = b„Q„ +l (x) + a„Q„(x) + c„Q„-i(x), n > 1 , (3.10) 

where b„, c, / 0 and c n /b n _\ > 0. Along with Q-i(x) — 0 and Qq(x) = 1, the 
three-term recurrence defines the polynomial system completely. Often the relation 
is written in a different form, 

Qn+i(x) = (A„x + B„)Q„(X) - C n Qn—i (^)j n > 0, 
and Favard proved the following converse result ([19]). 

Theorem 3.1 (Favard’s Theorem). Let A „ , B n , and C n be arbitrary sequences of 
real numbers and let { Q n (x)} be defined by the recurrence relation 

Qn+t (x) = (A„x + B n )Q„(x) - C n Q n - i(x), n > 0, 

together with Q o(x) = 1 and Q-\(x) = 0. Then the {Q n (x)} are a system of 
orthogonal polynomials if and only if A n ^ 0, C„ 0, and C„A„A n - 1 > 0 for 
all n. 
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3.1.3 Hypergeometric Series and the Askey Scheme 


Most orthogonal polynomials can be expressed in a unified way by using hyper- 
geometric series and incorporated in the Askey scheme. To this end, we first define 
the Pochhammer symbol (a) n as 

, _ f a, n — 0, 

{a)n ~ + 1)- ■■(« + «- 1), n = 1, 2, 

If a e N is an integer, then 

(a + n — 1)! 


(3.11) 


and for general a e 


(a)„ = 


C a)n = 


(a - 1)! 

T(a + n) 

r(fl) ’ 


n > 0, 


n > 0. 


The generalized hypergeometric series r F s is defined by 

(a\)k ■ ■ ■ ( a r )k z k 

{bi) k ---(b s ) k kl 


r F s (a \, . . . , a r \ b s ; z ) = 


(3.12) 


where bj ^0,-1, —2, . . . , for all i = 1, . . . , s. There are r parameters in the nu- 
merator and ,v parameters in the denominator. Clearly, the orders of the parameters 
are immaterial. 

For example, 


>« ; ^i = Es 


is the power series for an exponential function. 

When the series is infinite, its radius of convergence p is 


P = 


r < s + 1, 
r = s + 1, 
r > s + 1. 


If one of the numerator parameters a t , i = 1 r, is a negative integer, say, 

fli = —n, the series terminates because ( a\) k = (— «)t = 0,for£ = n + l,n-|- 
2, , and becomes 


*F,(a i 


,a r \b\, . . . , b s ; z) = 


( ~n)k ■ ■ ■ ( a r )k z k 

(/-iU •••'/.,)* k\ 


(3.13) 


This is a polynomial of degree n . 

The Askey scheme, which can be represented as a tree structure as shown in 
figure 3.1, classifies the hypergeometric orthogonal polynomials and indicates the 
limit relations between them. The tree starts with Wilson polynomials and Racah 
polynomials at the top. They both belong to class 4 Fj of hypergeometric orthogonal 
polynomials. Wilson polynomials are continuous polynomials and Racah polyno- 
mials are discrete. The lines connecting different polynomials denote the limit tran- 
sition relationships between them, which implies that polynomials at the lower ends 
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Figure 3. 1 The Askey scheme for hypergeometric orthogonal polynomials. 


of the lines can be obtained by taking the limit of some parameters in the polyno- 
mials at the upper ends. For example, the limit relation between Jacobi polynomials 
(x) and Hermite polynomials H n (x) is 


lim 

a— too 



2 n n\ ’ 


and that between Meixner polynomials M n (x; ft, c) and Charlier polynomials 
C n (x; a) is 


lim M„ ( x; B, ) = C n (x\ a). 

ft^oc V a + /3j 

For a detailed account of hypergeometric polynomials and the Askey scheme, the 
interested reader should consult [60] and [91]. 


3.1.4 Examples of Orthogonal Polynomials 

Here we present several orthogonal polynomials that will be used extensively in 
this book. The focus is on continuous polynomials. More specifically, we discuss 
Legendre polynomials defined on [—1, 1], Hermite polynomials defined on R, and 
Laguerre polynomials defined on [0, oo). These correspond to polynomials with 
support on a bounded interval (with proper scaling), the entire real line, and the 
half real line, respectively. 
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3. 1.4.1 Legendre Polynomials 
Legendre polynomials 


P n (x) = jF\ ( —n, n + 1; 1; 


1 — x 


satisfy 


and 


2n + 1 n 

P n+ 1 = xP„(x) P„-i(x), n > 0, 

n + 1 n + 1 




P n (x) P m (x)dx = 


2n + 1 


(3.14) 


(3.15) 


(3.16) 


Obviously the weight function in the orthogonality relation is a constant; 
i.e., w(x) = 1. The first few Legendre polynomials are 

Pq{x)=\, P\(x) = x, P 2 (x) = ^x 2 -]^ 

Legendre polynomials are a special case of Jacobi polynomials P„ a '^\x) with 
parameters a = ft = 0. The details for Jacobi polynomials can be found in 
appendix A. 


3. 1.4.2 Hermite Polynomials 


Hermite polynomials 


satisfy 



and 


H n+ i(x) = xH n (x ) — «//„_ i(x), n > 0, 



H m (x)H n (x)w(x)dx = n\S mn , 


where 


w(x) = 



The first few Hermite polynomials are 

Hq (x ) = 1 , //i(x) = x, H 2 (x) = x 2 — 1, H 3 (x) = x 3 - 3x, 


(3.17) 


(3.18) 

(3.19) 


Note that the definition of H n {x) here is slightly different from the classical one 
used in the literature. Classical Hermite polynomials H n {x) are often defined by 

H n+ i(x) = 2xH n (x) — 2n H n _\(x), 


n > 0, 
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and 



H n (x)H m (x)w(x)dx = 2 "n\8 mn . 


where w(x) = -fie x L 1 . The two expressions are off by a scaling factor. We em- 
ploy H„(x) here to facilitate the discussions associated with probability theory. 


3.1.43 Laguerre Polynomials 
Laguerre polynomials 

L { fi{x) = (a + 1)n i/q (-«; a + 1; x ) , a > -1, (3.20) 

w! 

satisfy 

(n + lli^^x) = (-x + 2n+a+l)i ) ( “ ) (x)-(«+a)T‘“_\(x), n > 0, (3.21) 

and 

[ L ( fi(x)L { “\x)w{x)dx = F(w +a + ^ (3.22) 

Jo «! 

where 

w(x) = e~ x x a . 

Note that the sign of the leading coefficients of Laguerre polynomials flips with 
increasing degree. 


3.2 FUNDAMENTAL RESULTS OF POLYNOMIAL APPROXIMATION 

Let P„ be the linear space of polynomials of degree at most n \ i.e., 

P„ = spanfx 4 : k = 0, 1, . . . , n). (3.23) 

We begin with a classical theorem by Weierstrass in approximation theory. 

Theorem 3.2 (Weierstrass). Let I be a bounded interval and let f e C°( 7). Then, 
for any e > 0, we can find n e N and p e P„ such that 

I f(x)~ p(x)\ < e, Vx e 1. 

We skip the proof here. Interested readers can find the proof in various books on 
approximation theory, for example, [18, 105, 106]. This celebrated theorem states 
that any continuous function in a bounded closed interval can be uniformly approx- 
imated by polynomials. From this theorem, a large variety of sophisticated results 
have emerged. A natural question to ask is whether, among all the polynomials of 
degree less than or equal to a fixed integer n, it is possible to find one that best 
approximates a given continuous function / uniformly in /. In other words, we 
would like to study the existence of <p n (f ) e P„ such that 
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This problem admits a unique solution, though the proof is very involved. An ex- 
tensive and general treatise on this subject can be found in [105]. The nth-degree 
polynomial <f>„ (f) is called the polynomial of best uniform approximation of / in 
7. Following theorem 3.2, one immediately obtains 

lim ||/ — </>«(/) lloo = 0. 
n— too 


3.3 POLYNOMIAL PROJECTION 


From here on we do not restrict ourselves to bounded intervals and consider the 
general cases of 7, where I = [— 1, 1], I = [0, oo[, or I — R. 

Another best approximation problem can be formulated in terms of norms other 
than the infinity norm used in (3.24). To this end, we define, for a positive weight 
function w(x), x e 7, the weighted L 2 space by 


|/ v 2 (x)w(x)dx < oo | (3.25) 


with the inner product 

M v) L 2 (/) = l u(x)v(x)vj(x)dx, Vw, v e 7,^(7 ), (3.26) 


and the norm 


Ml Ll(i) = (J u 2 {x)w(x)dx\ . (3.27) 


Throughout this book, we will often use the simplified notation ( u , v) w and \\u 
to stand for (m, v )^ ( /) and respectively, unless confusion would arise. 


3.3.1 Orthogonal Projection 

Let A be a fixed nonnegative integer and let {cpk(x)}^ =0 C P*/ be orthogonal poly- 
nomials of degree at most N with respect to the positive weight in (x ) ; i.e., 

(4>m(x), <t>n(x))Ll(I) = 110m IliJ (I) S m,n, 0 <m,tl <N. (3.28) 

We introduce the projection operator Pn-L 2 w (7) — > PA, such that, for any function 
/ e L 2 W (I), 

N 

PNf=J2fkMx), (3.29) 


fk = -V(/, 4>k)ii, o <k<N. (3.30) 

Obviously, P N f e P^. It is called the orthogonal projection of / onto P^ via 
the inner product (-,-)ij, and {/-} are the (generalized) Fourier coefficients. 
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The following trivial facts hold: 

P N f = I V/ e Pjv, 

P\i(pk = 0, Vk > N. 

Moreover, we have the following theorem. 

Theorem 3.3. For any f e i^(7) and any N e No, Pjv / best approximation 
in the weighted Lr norm (3.27) in the sense that 

\\f ~P N f hi = inf Wf-nLl- (3.31) 

Proof. Any polynomial i/r e Pjy can be written in the form xfr = 'f2k=o c k4>k for 
some real coefficients c*, 0 < k < N . Minimizing \\f — i/f|| l\ is equivalent to 
minimizing \\f — \j/ \\ 2 L2 , whose derivatives are 

9 d / N N \ 

-7-\\f-n 2 L 2 = — [\\ff L 2 -2j2c k (f,MLi + J2 c ^ k ^- 

dc J ” clc J \ " k = o *=o 7 

= -2 (f<Pj)Li +2cy||0y||^2, 0 < j < N. 

By setting the derivatives to zero, the unique minimum is attained when Cj = fj, 
0 < j < N, where f) are the Fourier coefficients of / in (3.30). This completes 
the proof. 

The projection operator also takes the name orthogonal projector in the sense 
that the error / — Py f is orthogonal to the polynomial space IP V . 

Theorem 3.4. For any f e L 2 W (I) and N e No, 

J(f - P N f)<pwdx = (/ - P N f 4>) Ll = 0, Vcf>e Py. (3.32) 

Proof. Let f eV y and define G : R. — > R. by 

G(v)±\\f-P N f+v4>\\ 2 Li , veR. 

From theorem 3.3, v = 0 is a minimum of G. Therefore, 

G'(v) = 2 J(f-P N f)fwdx + 2vU\\ 2 Ll 
should satisfy G'(0) = 0. And (3.32) follows directly. 

From (3.32) we immediately obtain the Schwarz inequality 

\\PNf\\Ll<\\f\\Ll (3.33) 

and the Parseval identity 

OO 

k = 0 


(3.34) 
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3.3.2 Spectral Convergence 

The convergence of the orthogonal projection can be stated as follows. 

Theorem 3.5. For any f e L 2 W (I), 

lim II/— Av/lliJ, = 0- (3.35) 

N^-OO 

We skip the proof here. When / is bounded, the proof is straightforward; see, for 
example, [34]. When 1 is not bounded, the proof is more delicate and we refer 
readers to [22] for details. 

The rate of convergence depends on the regularity of / and the type of orthogonal 
polynomials {</>*-}. There is a large amount of literature devoted to this subject. As 
a demonstration we present here a result for Legendre polynomials. 

Define a weighted Sobolev space for k = 0, 1, 2, . . . , by 


//*(/) = ju : / 

equipped an inner product 


d m v o 

- — e L 2 (I), 0 < m < k 
dx m 


A ^ — \ ( d m u d m v 

= E(^- 5 ^ 


(3.36) 


(3.37) 


and a norm ||«||//4 = (m, m)^*. 

Let us consider the case of / = [—1,1] with weight function w (x ) = 1 and 
Legendre polynomials \P n (x)} (section 3.1.4). The orthogonal projection for any 
f{x)eL 2 w (I) is 


p N m = Y J fkPkix), f k = 


i 


ii aid 


-(/; Pk)n- 


(3.38) 


k = 0 Ll 

The following result holds. 

Theorem 3.6. For any f{x) e H^[— 1, 1], p > 0, there exists a constant C, inde- 
pendent of N, such that 

||/- iWlkt-ui < CN-P\\f\\ m _ lM . (3.39) 

Proof. Since the Legendre polynomials satisfy (see (A.25) in appendix A) 

Q\Pk\ = k k P k , 

where 


Q = 


dx 


(1 -x 2 ) 


dx 


= (1-x 2 ) 


df 

dx 2 


2x- 


dx 


and X k = —k(k + 1). We then have 


(/, Pk)il = ~ f Q[Pk]f(x)dx = 1 [ ((T - x 2 ) Pi'f - 2 xP'f) dx 
A'k J - 1 J - 1 

1 c l 

= / [((1 - X 2 )f ) ’Pi + 2* A/] dx, 

^k J - 1 
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where integration by parts has been applied to the first term of the integrand to 
derive the last equality. Upon simplifying the last expression, we obtain 

(/, P k ) Li = -- I (1 - x 2 )f'P'dx = 1 f ((1 - x 2 )/)' P k dx, 

J - 1 J - 1 

where integration by parts is again utilized. This implies 

(/; Pk),i = T(g[/],p ,) i2 
4 k 

By applying the procedure repeatedly for m times, we obtain 


(./; Pk)ii = T (Q"‘[f], p k ) L1 . 

k k 

The projection error can be estimated as 


\f-p N f\\ 2 Ll = E fk\\ p k\\li = E 


1 


k=N+l 
oo 

= E 

k=N + 1 ''’k 
< X 


II Pk II 


-ifPkY, 


II Pk 


k "Ll 


k=N + 1 

-iorin, Pk)% 


Ll 


N 


oo ^ 

E^(e m [/]^4 


k=0 


< N \\Q m [f]\\ L i i 
<CN~ Am || f\ |L„, 


where the last inequality relies on ||2 m [/]lliJ, — C||/|| // 2 m, which is a direct conse- 
quence of the definitions of Q m [f \ and the norms. By taking p — 2m, the theorem 
is established. 


Therefore, the rate of convergence of the Legendre approximation relies on the 
smoothness of the function /, measured by its differentiability. For a fixed approx- 
imation order N, the smoother the function /, the larger the value of p. and the 
smaller the approximation error. This kind of convergence rate is referred to in the 
literature as spectral convergence, ft is in contrast to the traditional finite difference 
or finite element approximations where the rate of convergence is fixed regardless 
of the smoothness of the function. An example of spectral convergence is shown in 
figure 3.2, where the error convergence of the Legendre projections of |sin(^x)| 3 
and |x| are given. Both functions have finite, but different, smoothness. The con- 
vergence rates of the two functions are clearly different in this log-log figure, with 
| sin(7 T jc) | 3 having a faster rate because of its higher differentiability. 

If f(x) is analytic, i.e., infinitely smooth, the convergence rate is faster than any 
algebraic order and we expect 

U-Pn/UlI ~Ce-“ JV ||/|| i 2, 

where C and a are generic positive constants. Thus, for an analytic function, spec- 
tral convergence becomes exponential convergence. An example of exponential 
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Figure 3.2 Spectral convergence: projection error of |sin(7rx)| 3 and |x| by Legendre poly- 
nomials in x e [—1, 1], (N is the order of expansion.) 



Figure 3.3 Exponential convergence: projection error of cos(rrx) by Legendre polynomials 
in x e [—1, !]• (N is the order of expansion.) 


convergence is shown in figure 3.3, where the projection error of cos(ttx) by Legen- 
dre polynomials is plotted. Exponential convergence is visible in this kind of semi- 
log plot. 


3.3.3 Gibbs Phenomenon 

When the function / is not analytic, the rate of convergence of the polynomial 
projection is no longer faster than the algebraic rate. In the case of discontinuous 
functions, the convergence rate deteriorates significantly. For example, consider the 
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Figure 3.4 Orthogonal series expansion of the sign function by Legendre polynomials. 
(N is the order of expansion.) 


sign function in (—1, 1): 


sgn(x) = 



x > 0; 
x < 0. 


The Legendre series expansion is 


OO 

sgn(x) = Y] 

n = 0 


(— 1)"(4« + 3) (2m)! 
2 2 ”+ 1 0 + 1)!«! 


^2«+l(x). 


The partial sums of the series are plotted in figure 3.4 for several values of N. 
We observe oscillations near the discontinuity, and they do not disappear as N is 
increased. This is referred to as the Gibbs phenomenon. It is a numerical artifact of 
using globally smooth polynomial basis functions to approximate a discontinuous 
function. In fact, for this Legendre series expansion, the Gibbs phenomenon has 
a long-range effect in the sense that it seriously affects the rate of convergence at 
the endpoints x = ± 1 of the interval. For more a detailed discussion of the Gibbs 
phenomenon, see [47], 


3.4 POLYNOMIAL INTERPOLATION 

The goal of polynomial interpolation is to construct a polynomial approximation to 
a function whose values are known at some discrete points. More precisely, given 
m + 1 pairs of (x,-, y, ) , the problem consists of funding a function G = G(x) such 
that G(xj) = yt for i = 0. . ... in, where y, are given values and G interpolates {>»,•} 
at the nodes {x,}. In polynomial interpolation, G can be an algebraic polynomial, 
a trigonometric polynomial, a piecewise polynomial (that is, a local polynomial), a 


POLYNOMIALS AND APPROXIMATION 


37 


rational polynomial, etc. In the brief introduction here, we focus on global polyno- 
mials of algebraic form. 

3.4.1 Existence 

Let us consider TV + 1 pairs of (x , , y, ) . The problem is to find a polynomial Q ;W e 
Pm, called an interpolating polynomial, such that 

Q M (xi) — ctMxf 1 H + a\Xi + ao = i = 0, . . . , TV. (3.40) 

The points {x,} are interpolation nodes. If TV ^ A/, the problem is over- or under- 
determined. If TV = M , the following results hold. 

Theorem 3.7. Given TV + 1 distinct points xo, . . . , x.y and TV + I corresponding 
values _vo, . . . , yxt, there exists a unique polynomial Q ^ e Py such that Q^(xj) = 
Vi for i = 0, . . . , TV. 

Proof. To prove existence, let us use a constructive approach that provides an ex- 
pression for Qn. Denoting {/, }^ 0 as a basis for Pjy, then Qn(x) = ^ /=0 bjlj(x) 
with the property that 

N 

Qn{x{) = YbjljlXj) = Vi, 

7=0 

Let us define 

N+\ 

l, e P,v : /,(,)= 

1 L Xi — X 
7=0 
j& 

then l,(xj) = Sij and we obtain bi = yt. 

It is easy to verify that {/,, i — 0, .... TV} form a basis for Pjy (left as an ex- 
ercise). Consequently, the interpolating polynomial exists and has the following 
form, called the Lagrange form , 

N 

Q N (.x) = Y j yi l i(.x). (3.42) 

i=0 

To prove uniqueness, suppose that another interpolating polynomial Qm(x) of 
degree M < N exists such that Qm(x ,) = y;, i = 0, . . . , N. Then, the difference 
polynomial Q N — Q M vanishes at TV + 1 distinct points x, and thus coincides with 
the null polynomial. Therefore, Qm — Qn- 

Another approach also provides a way of constructing the interpolating polyno- 
mial. Let 

Q N (x.i) — arjx* 1 H + a\x t + a 0 = J7, i = 0, . . . , N. (3.43) 

This is a system of TV + 1 equations with TV + 1 unknowns of the coefficients 
ao,---, a N . By letting a = (a 0 , ■■■, a N ) T , y = (y 0 , . . . , y N ) T , and A = (a, 7 ) = 
(x/), system (3.43) can be written as 


i = 0, . . . , TV. 


i = 0 TV; (3.41) 


Aa = y. 
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Figure 3.5 Polynomial interpolation of /(x) = 1/(1 +25x 2 ) on [— 1, 1], the rescaled Runge 
function. Left: interpolation on uniformly distributed nodes. Right: interpolation 
on nonuniform nodes (the zeros of Chebyshev polynomials). 


The matrix A is called a Vandermonde matrix, which can be shown to be nonsin- 
gular (left as an exercise). Therefore, a unique set of solutions to the coefficients a 
exists. 


3.4.2 Interpolation Error 

Let f(x),x e I, be a given function and let LLy/ ( x ) be its interpolating polynomial 
of the ATh degree constructed by using the values of f (x) at N + 1 distinct points. 
Then the following result holds. 

Theorem 3.8. Let xo , . . . , xjy be N + 1 distinct nodes and let x be a point inside the 
interval I. Assume that f G C N+l (I x ), where I x is the smallest interval containing 
the nodes xo, . . . , Xjv and x. Then the interpolation error at the point x is 

E n (x) = fix) - n N f(x ) = (3-44) 

where § G I x and q^+\ = flS'C* — x i) ‘ s the nodal polynomial of degree N + 1. 

The proof for this standard result is skipped here. (See, for example, [6]). 

We should note that a high-degree polynomial interpolation with a set of uni- 
formly distributed nodes is likely to lead to problems, with qu+\ ix) behaving rather 
wildly near the endpoint nodes. This leads to ITjv fix) failing to converge for sim- 
ple functions such as fix) = (1 + x 2 ) 1 on [—5, 5], a famous example due to Carl 
Runge termed the Runge phenomenon. To circumvent the difficulty, it is essential to 
utilize piecewise low-degree polynomial interpolation or high-degree interpolation 
on nonuniform nodes. In the latter approach, the zeros of orthogonal polynomi- 
als provide an excellent choice. This can be seen in figure 3.5, where the Runge 
function is rescaled to domain [—1, 1] and interpolated by 26 points (A = 25). 
The interpolation on uniformly distributed nodes (on the left in figure 3.5) becomes 
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ill-conditioned and incurs large errors. The interpolation on nonuniformly distri- 
buted nodes (on the right in figure 3.5) is stable and accurate. 


3.5 ZEROS OF ORTHOGONAL POLYNOMIALS AND QUADRATURE 

It is well known that a polynomial of degree N has at most N distinct complex ze- 
ros. Although very little can be said about real zeros in general cases, the following 
result holds. 

Theorem 3.9. Let { Q n (x)}„ g u, x e 1, be orthogonal polynomials satisfying ortho- 
gonal relation (3.8). Then, for any n > 1, Q„ has exactly n real distinct zeros in I. 

Proof. First note that ( Q„ , 1 )„ : = 0. Therefore, Q„ changes sign in 7, hence it has 
at least one real zero z\ e I. And n — I is proved. For n > 1, we can find another 
zero zi £ / with zi f z\ since if Q„ vanishes only at z\, then the polynomial 
(x — z\ )Q n would not change sign in 7, which is in contradiction to the relation 
((x — z\), Q„) w = 0, obtained by orthogonality. In a similar fashion, we consider 
the polynomial (x — z\)(x — zf)Q n and, if n >2, deduce the existence of a third 
zero, and so on. The procedure ends when all n zeroes are obtained. 

There are many properties of the zeros of orthogonal polynomials. For example, 
one can prove that Q n and Q n -\ do not have common zeros. Moreover, between 
any two neighboring zeros of Q n -\, there exists one and only one zero of Q n . Let 
us recall the following statement. 

Theorem 3.10. Let { Q„ }„ g H be a sequence of orthogonal polynomials in I. Then, 
for any interval [a, b] C I, a < b, it is possible to find m e N such that Q m has at 
least one zero in [a, b]. 

In other words, this theorem states that the set J = U„>i Ua-=i{ z [" ) } is dense in 
7, where {z < / " > } are the zeros of the orthogonal polynomials Q n . The proof can be 
found in classical texts on polynomials such as [100]. 

For example, the zeros of the Legendre polynomials P n (x), x £ [-1, 1], satisfy 

k — \ k 

— 1 < — cos r < < — cos T zt <1, 1 < k < n. 

n+\ n+\ 

The length of the interval between two consecutive zeros is 

k k — i 2k — j \ 

L = — cos T n + cos T7r = 2 sin -it sin — - — tc. 

n + \ n+\ 2n + l 2n + 1 

For large n, L oc n~ 2 for k ^ 1 or k ~ n; and L oc n~ l for moderate values 
of k. This indicates that the zeros of Legendre polynomials are clustered toward the 
endpoints of the interval [—1, 1], a feature shared by many orthogonal polynomials. 

The nonuniform distribution of the zeros of orthogonal polynomials makes them 
excellent candidates for polynomial interpolation (see figure 3.5). In addition, they 
are also excellent for numerical integration. 
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Let 

Ilf] = J f(x)u>(x)dx (3.45) 

and define an integration formula with q > 1 points, 

q 

W[f] = J2 f(x (j) )w U) , (3.46) 

j = i 

where x (j> are a set of nodes and w (j> are integration weights, j = 1, . . . , q. The 
objective is to find a set of {x (y ) ; w^} such that U q [f] /[/] and hopefully 
lim^^oo U q [ f] = /[/]. For example, the well-known trapezoidal rule approxi- 
mates (3.45) in an interval [a, b ] by 

nf]^ b -^ L im + /m, 

which implies that, in (3.46), q = 2, {x (y, J = {a, b}, and = ff-, j = 1, 2. 

Highly accurate integration formulas can be constructed by using orthogonal 
polynomials. Let {(2«l«eN be orthogonal polynomials satisfying (3.8) and let 
{ z A- W, }/tLi b e the zeros of Qn- Let /* :V 11 be the (N — l)th-degree Lagrange polyno- 
mials through the nodes 1 < j < N, and let FIjv-i/(x) = ^^ =1 f(z { J N) )l { J N ~ l) 
be the (N — l)th-degree interpolation of f(x). Then, the integral (7.12) can be ap- 
proximated by integrating n v _ i f(x), 

r N 

/ f(x)w(x)dx & y] f(zf ] )w\ N \ (3.47) 

Jl 7=1 

where 

w ( P = Jlf^wdx, 1 <j<N, 

are the weights. This approximation is obviously exact if / e Pjv-i- However, the 
following result indicates that it is more accurate than this. 

Theorem 3.11. Formula (3.47) is exact; i. e. , it becomes an equality if f(x) is any 
polynomial of degree less than or equal to 2N — 1 in I. 


Proof For any / e P 2 W- 1 , let q = Yl^-if e Pv-i be the (N — l)th-degree 
interpolation using {z ( f Then /— q vanishes at and can be expressed 

as / — q = Q N (x)r{x), where r(x) is a polynomial of degree at most N — I . By 
orthogonality and using the exactness of (3.47) for polynomials of degree up to 
N — 1 , we have 


J fwdx = J qwdx + J Q^rwdx = J qwdx 

=x:,(zf ) ^=x:/(zf)u ) f. 

7=1 7=1 
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The converse of the theorem is also true; i.e., if formula (3.47) holds for any 
/ e IP 2 /V - 1 ■ then the nodes are zeros of (Ty The degree of the integration formula 
cannot be improved further. In fact, if one takes / = Q 2 N e IVy, the right-hand 
side of (3.47) vanishes because Qn(z\ N> ) = 0 , j = 1, . . . , N, but the left-hand 
side does not. 


3.6 DISCRETE PROJECTION 


Here we define a discrete projection of a given function / e L„,(/) as 

N 

w(x) = X!uw. (3- 48 ) 

n = 0 

where the expansion coefficients 

~ 1 1 J 1 

fn = — £/ ? [/(x)0„(x)] = — Y f(x u) )(p n (x {j) )w i, 0 < n < N. 

\\M 2 Ll WMliU ~ ~ 

(3.49) 

When an integration formula JJ q is used, the coefficients f n are approximations of 
the coefficients f n (3.30) in the continuous orthogonal projection (3.29). That is, 

f n Kif n =J f(f) n wdx, 0 < n < N . 

Note that this definition is slightly more general than what is often used in the 
literature, in the sense that we did not specify the type and number of nodes used 
in the integration formula. The only requirement is that the integration formula 
approximates the corresponding continuous integrals. 

In classical spectral methods analysis, the integration formula is typically chosen 
to be the Gauss quadrature corresponding to the orthogonal polynomials {</>„(x)J 
satisfying (3.28). More specifically, let {z^ V| } ^| 1 be the zeros of </>v+ i(x) in I. Let 
Thv/be the Lagrange interpolation of / (x) in the form of (3.42), 

N + 1 

n N f(x) = y nz^ijix). (3.50) 

/—i 

Let us now use the (N + l)-point Gauss quadrature to evaluate the coefficients in 
the discrete expansion (3.48). That is, 


fn = 


N+ 1 


.Wk „(AD 


0 <11 <N. (3.51) 


Then, the following results hold. 

Theorem 3.12. Let 1^ f be defined by (3.48), where the coejficients {f n } are eval- 
uated by ( N + \)-point Gauss quadrature based on orthogonal polynomial {</>„}, 
as in (3.51). Let IT^ f be the Lagrange interpolation of f through the same set of 
Gauss nodes, as in (3.50). Then, for any f e Pv, TIjv f = In /■ 
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Proof. For / E Pj, n v f = / because of the uniqueness of Lagrange interpola- 
tion. The coefficients of its discrete expansion are 


, AM-1 




Pn(z )W 


*•"11 7=1 


ll&ll?2 l 


= fn , 


f(x)<f> n (x)w(x)dx 
0 < n < N , 


because the (N + l)-point Gauss quadrature is exact for integrating polynomials 
of degree up to 2N + 1 and the integrand is in P 2 N. Therefore, hf = P N f, 
the continuous orthogonal projection of /, which in turn equals /. The proof is 
established. 


The above equivalence of the discrete projection, continuous projection, and La- 
grange interpolation does not hold for a general function / because the Gauss 
quadrature will not be exact. In such a case, f n ^ f n . Hence, Pm f f In /■ The dif- 
ference between the continuous orthogonal projection and the discrete projection is 
often termed aliasing error, 

A N f=P N f-I N f (3.52) 

By letting y„ = || || and realizing that the summation in (3.51) defines a 
discrete inner product [/, <p n ] w , we obtain, for all 0 < n < N, 


fn = ~[f Mu, = ~ 

Yn Vn 




j = 0 


j OO 

= ]» 


y=o 


1 / °° 

- — E^’ 

Yn \]<N j>N 

J OO 


n\w fj 


j>N 


Therefore, 


N 

A N f=J2^n-fn)4h, 

n= 0 

N ^ 00 00 N ^ 

= E 7 E fMp = E E —Mb’ Mwfj 

n= 0 ]>N j>N n= 0 


j>N 

Thus, the aliasing error can be seen as the error introduced by using the interpola- 
tion of the basis, I]f<f>j, rather than the basis itself to represent the higher expansion 
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modes (j > N ). The aliasing error stems from the fact that one cannot distinguish 
between lower and higher basis modes on a finite grid. A general result holds that 
the aliasing error induced by Gauss points is usually of the same order as that of the 
projection error. Hence for well-resolved smooth functions, the qualitative behavior 
of the continuous and the discrete expansions is similar for all practical purposes. 
We will not pursue a further in-depth discussion of this and instead refer interested 
readers to the literature, for example, [51]. 


Chapter Four 


Formulation of Stochastic Systems 


This chapter is devoted to the general aspects of formulating stochastic equations, 
i.e., given an established deterministic model for a physical system, how to prop- 
erly set up a stochastic model to study the effect of uncertainty in the inputs to the 
system. Prior to any simulation, the key step is to properly characterize the ran- 
dom inputs. More specifically, the goal is to reduce the infinite-dimensional prob- 
ability space to a finite-dimensional space that is amenable to computing. This is 
accomplished by parameterizing the probability space by a set of a finite num- 
ber of random variables. More importantly, it is desirable to require the set of 
random variables to be mutually independent. We remark that the independence 
requirement is very much a concern from a practical point of view, for most, if 
not all, available numerical techniques require independence. This is not as strong 
a requirement from a theoretical point of view. Although there exist some tech- 
niques to loosen it, in this book we shall continue to employ this widely adopted 
requirement. 

To summarize, the critical step in formulating a stochastic system is to properly 
characterize the probability space defined by the random inputs by a set of a finite 
number of mutually independent random variables. In many cases such a charac- 
terization procedure cannot be done exactly and will induce approximation errors. 


4.1 INPUT PARAMETERIZATION: RANDOM PARAMETERS 

When the random inputs to a system are the system parameters, the parameteriza- 
tion procedure is straightforward, for the inputs are already in the form of parame- 
ters. The more important issue is then to identify the independent parameters in the 
set. The problem can be stated as follows. 

Let Y = (Yi , . . . , Y„), n > 1, be the system parameters with a pre- 
scribed distribution function Fy (_y) = P(Y < y), y e R' ! , and find 
a set of mutually independent random variables Z = (Z\. . . . , Zf) e 
W 1 , where 1 < d < n, such that Y = T(Z) for a suitable transforma- 
tion function T . 

Let us use a simple example to illustrate the idea. Consider an ordinary differen- 
tial equation with two random parameters, 


— (t, co) = — a(co)u , 
dt 


u( 0, co) = /3(<y), 


(4.1) 
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where the rate constant a and the initial condition /I are assumed to be random. 
Thus, the input random variables are Y ( a > ) = (a, ft) e R 2 . 

If a and ft are mutually independent, then we simply let Z(a>) = Y(co). The 
solution 


u(t, a>) : [0, T] x £2 — »■ R. 

can now be expressed as 

u{t, Z) : [0, f]xI 2 -> R, 

which has one time dimension and two random dimensions. 

If a and ft are not independent of each other, it implies that there exists a function 
/ such that 

/(«, ft) = 0 . 

Then it is possible to find a random variable Z (o>) to parameterize the relation such 
that 


a(a > ) = a(Z(co)), ft(a>) — b(Z(a>)), 

and f{a, b) = 0. Or, equivalently, the dependence between a and ft implies that 
there exists a function g such that 

ft = g( a). 

Then we can let Z(co) = a(o>) and ftioj) = g(Z(co)). Which approach is more 
convenient in practice is problem-dependent. Nevertheless, the random inputs via 
a and ft can now be expressed via a single random variable Z, and the solution 
becomes 


u{t, Z) : [0, T] x R — ► R, 

which now has only one random dimension. 

In practice, when there are many parameters in a system, finding the exact form 
of the functional dependence among all the parameters can be a challenging (and 
unnecessary) task. This is especially true when the only available information is the 
(joint) probability distributions of the parameters. The goal now is to transform the 
parameters to a set of independent random parameters by using their distribution 
functions. 

4.1.1 Gaussian Parameters 

Since the first two moments, mean and covariance, can completely characterize 
Gaussian distribution, the parameterization problem can be solved in a straightfor- 
ward manner. 

Let Y = (Ti, . . . , Y„) be a random vector with a Gaussian distribution of 
J\f(0, C), where C e R" x,! is the covariance matrix and the expectation is assumed 
to be zero (without loss of generality). Let Z ~ A r (0, 1), where I is the n x n 
identity matrix, be a uncorrelated Gaussian vector of size n. Thus, the components 
of Z are mutually independent. Let A be an n x n matrix; then by theorem 2.16, 
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A Z ~ Af (0, AA 7 ). Therefore, if one finds a matrix A such that AA 7 = C, then 
Y — A Z will have the given distribution Af( 0, C). This result is a special case of a 
more general theorem by Anderson [4] . 

Since C is real and symmetric, solving the problem AA r = C can be readily 
done via, for example, Cholesky’s decomposition, where A takes the form of a 
lower-triangular matrix: 


an = cn/JcTu !</<«, 

an = yj cu - Y!k=\ af k , 1 <i<n, 

a ij = (cij - i a ‘k a jk) /ajj, 1 < j < i <n, 


(4.2) 


di j — 0 , 


i < j < n , 


where ajj and c ly , 1 < i, j < n, are the entries for the matrices A and C, respec- 
tively. 


4.1.2 Non-Gaussian Parameters 

When the system parameters have a non-Gaussian distribution, the parameteri- 
zation problem is distinctly more difficult. However, there exists a remarkably 
simple transformation due to Rosenblatt [88] that can accomplish the goal. Let 
Y = (Ti, . . . , Y„) be a random vector with a (non-Gaussian) distribution function 
F Y (v ) = P(Y < y) and let z = (z 1; . . . , z„) = Ty = T (yi, . . . , y n ) be a transfor- 
mation defined as 

zi = P(Y 1 <y l ) = F l (y l ), 

Z2 = P(Y 2 < yi | Y\ = y\) = F 2 (yi\y\), 

(4.3) 


Z* = P(Yn < Yn I Y n —i = V„- 1 , . . . , Ti = jV) = F„(y n \y n -u . . . , Vi). 

It can then be shown that 

P(Zi < z,- i = 1 ,...,«) 

J d y „ F n (y„\y„-i, . . . , y\) ■ ■ ■ d yi F, (y t ) 

* Zl 11 

dz\ ■ ■ ■ dz n = J~ | z,-, 

1 i=i 

where 0 < z,- < 1 , i = 1 Hence Z = (Z \, . . . , Z n ) are independent and 
identically distributed (i.i.d.) random variables with uniform distribution in [0, 1]". 

Though mathematically simple and powerful, Rosenblatt transformation is not 
easy to carry out in practice, for it relies on the conditional probability distributions 
among the random parameters. Such information is rarely known completely. And 
even if it is known, it is rarely given in explicit formulas. In practice, some kinds of 
numerical approximations of Rosenblatt transformation are required. This remains 
an understudied topic and is beyond the scope of this book. 
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4.2 INPUT PARAMETERIZATION: RANDOM PROCESSES AND 
DIMENSION REDUCTION 

In many cases, the random inputs are stochastic processes. For example, the inputs 
could be a time-dependent random forcing term that is a stochastic process in time, 
or an uncertain material property, e.g., conductivity, that is a stochastic process in 
space. The parameterization problem can be stated as follows. 

Let (Y t , t e T) be a stochastic process that models the random in- 
puts, where 1 is the index belonging to an index set T , and find a 
suitable transformation function R such that Y, = R(Z), where Z = 

(Zi, . . . , Zj), d > 1, are mutually independent. 

Note that the index set T can be in either the space domain or the time domain 
and is usually an infinite-dimensional object. Since we require d to be a finite inte- 
ger, the transformation cannot be exact. Therefore, Y, ~ R(Z) in a proper norm or 
metric, and the accuracy of the approximation will be problem-dependent. 

A straightforward approach is to consider the finite-dimensional version of Y, 
instead of Y, directly. This requires one to first discretize the index domain T into 
a set of finite indices and then study the process 


(Y, YJ, ti, ... ,t„ e T, 


which is now a finite-dimensional random vector. The parameterization techniques 
for random parameters from the previous section can now be readily applied, e.g., 
Rosenblatt transformation. 

The discretization of the Y, into its finite-dimensional version is obviously an 
approximation. The finer the discretization, the better the approximation. However, 
this is not desired because a finer discretization leads to a larger dimension of n 
and can significantly increase the computational burden. Some kinds of dimension 
reduction techniques are required to keep the dimension as low as possible while 
maintaining a satisfactory approximation accuracy. 

4.2.1 Karhunen-Loeve Expansion 

The Karhunen-Loeve (KL) expansion (cf. [73], for example) is one of the most 
widely used techniques for dimension reduction in representing random processes. 

Let iiy(t) be the mean of the input process Y, and let C(t, s ) = cov(7 r , 7 S ) be 
its covariance function. The Karhunen-Loeve expansion of 7, is 


OO 



(4.4) 


where i /r,- are the orthogonal eigenfunctions and /., are the corresponding eigenval- 
ues of the eigenvalue problem 

J C(t, s)fi(s)ds = Xifiit), t e T, (4.5) 

and {7/(<y)} are mutually uncorrelated random variables satisfying 
E[7,] = 0, E [Y i Yj] = 8 ij , 


(4.6) 
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and defined by 

Yi(co) = -j= f T (Y t (fi>) - Vi. (4.7) 

The Karhunen-Loeve expansion, in the form of the equality (4.4), is of little use 
because it is an infinite series. In practice, one adopts a finite series expansion, e.g., 

d 

Y t {oo ) ~ Aiy(f) + E y/kiMOYtico), d> 1. (4.8) 

i= 1 

The natural question to ask is when to truncate the series. That is, how to choose 
d so that the approximation accuracy is satisfactory. The answer to the question is 
closely related to an important property of the Karhunen-Loeve expansion — decay 
of the eigenvalues A.,- as index i increases. Here we illustrate the property with the 
following examples. 


Example 4.1 (Exponential covariance function). Let C(t, s ) = exp(— \t — s\/a), 
where a > 0 is the correlation length, and let t e T = \—b. b\ be in a bounded do- 
main with length 2b. Then the eigenvalue problem (4.5) can be solved analytically 
[109]. The eigenvalues are 


ki = 


2 a 

1 +a 2 wf ’ 


2 a 

l+a 2 u? ’ 


if i is even, 
if i is odd, 


(4.9) 


and the corresponding eigenfunctions are 

if i is even, 
if i is odd, 

where w,- and v, are the solutions of the transcendental equations 

I aw + tan(u>6) = 0, for even i, 

1 — av tan(uf)) = 0, for odd i. 


sin(Wit) j singiM 
cos iVit) /Jb+^, 


(4.10) 


In figure 4. 1 , the first four eigenfunctions are shown for the exponential covariance 
function in [—1, 1]. It is obvious that the higher modes (the eigenfunctions with 
larger index i ) have a finer structure compared to the lower modes. The eigenvalues 
are shown in figure 4.2 for several different correlation lengths a. It can be seen 
that the eigenvalues decay, and the decay rate is larger when the correlation length 
is longer. When the correlation length is very small, e.g., a = 0.01, the decay of 
the eigenvalues is barely visible. 


Example 4.2 (Uncorrelated process). The limit of diminishing correlation length 
is the zero correlation case, when the covariance function takes the form of a delta 
function, C(t, s) = S(t — s). It is easy to see from (4.5) that now any orthogonal 
functions can be the eigenfunctions and the eigenvalues are a constant, i.e., = 1 , 

Vi. In this case, there will not be any decay of the eigenvalues. 
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Figure 4. 1 The first four eigenfunctions of the exponential covariance function. 



Figure 4.2 The first 20 eigenvalues of the exponential covariance function with different 
correlation lengths a. 


Example 4.3 (Fully correlated process). The other limit is when C(t, s ) = 1, 
which implies an infinite correlation length and that the process Y, is fully corre- 
lated. This is the rather trivial case where the process depends on just one random 
variable. It is straightforward from (4.5) to show that there exists one nonzero eigen- 
value corresponding to a constant eigenfunction and that the rest of the eigenvalues 


are zero. 
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The aforementioned examples illustrate an important property of the Karhunen- 
Loeve expansion: for a given covariance function, the decay rate of the eigenvalues 
depends inversely on the correlation length. Long correlation length implies that 
the process is strongly correlated and results in a fast decay of the eigenvalues. The 
limit of this, infinitely long correlation length, is the fully correlated case where the 
eigenvalues decay to zero immediately. Conversely, a weakly correlated process has 
short correlation length and results in a slow decay of the eigenvalues. The limit of 
this, the uncorrelated process with zero correlation length, has no eigenvalue decay. 

The decay rate of the eigenvalues provides a guideline for truncating the infinite 
KL series (4.4) into the finite KL series (4.8). The common approach for truncation 
is to examine the decay of A.,- and keep the first d eigenvalues so that the contribu- 
tion of the rest of the eigenvalues is negligible. How much is considered negligible, 
usually given in terms of a small percentage as a cutoff criterion, is specified on 
a problem-dependent basis. Naturally, for a given cutoff criterion, a strongly cor- 
related process allows a finite KL expansion with a fewer number of terms than a 
weakly correlated process does. 

There are many more properties regarding the KL expansion. For example, the 
error of a finite-term expansion is optimal in terms of the mean-square error. We 
will not devote further discussion to this and refer interested readers to references 
such as [93]. 


4.2.2 Gaussian Processes 

The truncated Karhunen-Loeve series (4.8) provides a way to approximate a ran- 
dom process by a function (a series) involving a finite number of random variables. 
The set of random variables Y, (o>) are uncorrelated, as in (4.6). For Gaussian ran- 
dom variables, uncorrelation and independence are equivalent. Furthermore, linear 
combinations of Gaussian random variables remain Gaussian-distributed. There- 
fore, if Y, (cu) is a Gaussian process, then the random variables Y, in (4.4) and (4.8) 
are independent Gaussian random variables. Hence (4.8) provides a natural way 
to parameterize a Gaussian process by a finite number of independent Gaussian 
random variables. 


4.2.3 Non-Gaussian Processes 

When the input random processes are non-Gaussian, their parameterization and di- 
mension reduction are significantly more challenging. The main problem is that, for 
non-Gaussian distributions, uncorrelation of the random variables Y, in (4.4) does 
not imply independence. Hence the Karhunen-Loeve expansion does not provide 
a way of parameterization with independent variables. In many practical computa- 
tions, one often still uses the KL expansion for an input process and then further 
assumes that the Y, are independent. Though this is not a rigorous approach, at the 
moment there are not many practical methods for achieving the parameterization 
procedure. We will not engage in further discussion of this because it remains an 
active (and open) research topic. 
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4.3 FORMULATION OF STOCHASTIC SYSTEMS 

We now illustrate the main steps in formulating a stochastic system by taking into 
account random inputs to a well-established deterministic system. We choose par- 
tial differential equations (PDEs) as a basic model, although the concept and pro- 
cedure are not restricted to PDEs. 

Let us consider a system of PDEs defined in a spatial domain D C R\ t = 
1, 2, 3, and a time domain [0, T] with T > 0, 

' u, (x, t, co) = C(u), Dx(0,T]x£i, 

B(u) = 0, 35 x [0, f] x (4.11) 

u — uq, D x {t = 0} x £2, 

where £ is a (nonlinear) differential operator, B is the boundary condition operator, 
mo is the initial condition, and a> e £2 denotes the random inputs of the system 
in a properly defined probability space (f2, T , P ). Note in general that it is not 
important, nor is it relevant, to identify precisely the probability space. The solution 
is therefore a random quantity, 

u{x, t, cu) : D x [0, T] x £2 — > R"“, (4.12) 

where n u > 1 is the dimension of u. 

The random inputs to (4. 1 1) can take the form of random parameters and random 
processes. Let us assume that they can all be properly parameterized by a set of 
independent random variables using the techniques discussed in the previous two 
sections. Let Z = (Z\, . . . , Zj) e W l , d > 1, be the set of independent random 
variables characterizing the random inputs. We can then rewrite system (4.1 1) as 

' u,(x, t, Z) = C(u), D x (0, T] x 

B(u) = 0, dD x [0, T] x M 0 ', (4.13) 

U = Uq, D X {t — 0} X 

The solution is now 

u(x, t, Z) : D x [0, T]x IT' -* R"“. (4.14) 

The fundamental assumption we make is that (4.11) is a well-posed system P- 
almost surely in Q. Loosely and intuitively speaking, this means that if one gener- 
ates an ensemble of (4.13) by generating a collection of realizations of the random 
variables Z, then each realization is well posed in its corresponding deterministic 
sense. 

Example 4.4. Consider the same example of the ordinary differential equation 
(ODE) (4.1) 

du 

— (t, u > ) = —a{cd)u, «( 0, co) = P(co). 

dt 


(4.15) 
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If the input random variables a and ft are independent, then we let Z = (Z i , Zft) = 
(a, ft) and rewrite the problem as 
du 

— (t, aft) = —Z\u, u(0,co) = Zi. (4.16) 

dt 

If, however, a and ft are dependent, then, as discussed in section 4.1, we can 
let Z — a and there exists a function such that ft — g(a). The system can be 
rewritten as 

— ( t,co) = —Zu , w(0, co) = g(Z). (4.17) 


Example 4.5 (Stochastic diffusion equation). Consider a one-dimensional sto- 
chastic elliptic equation 


f V • (k(x, w)Vw) = f(x,co), x G (— 1, 1), 
[w(— 1, co) — ui(co), u(\, co) = u r (co), 


(4.18) 


where the diffusivity field k and the source term / are assumed to be random fields 
and Ui and are random variables. For simplicity of exposition, only the Dirichlet 
boundary condition is considered. Let us assume that the diffusivity field k can be 
parameterized by a truncated KL expansion (4.8) with d K terms, i.e., 


<4 

k(x, a>) ~ /c(x, Z K ) = /Xk(x) + ^^Kj(x)Zf (co), 


where the functions £, (x) are determined by the eigenvalues and eigenfunctions of 
the covariance function of k(x, oo), and Z) (co) are mutually independent. Similarly, 
let f(x,co) be parameterized as 


d f 

f(x, co) « f(x, Z f ) = fi f (x) + Y2 fj (x ) z/ (co), 


with d f terms and mutually independent zj (co) . Let us assume k and / are inde- 
pendent of each other and also independent of m and u r . Then let 

Z = (Z 1 , ...,Z d )= (Zf, ...,Z K dK ,z{,..., z{ f , Ui, Ur), 
where d = d K + df + 2, and we can write the elliptic problem as 
[V-(k(x,Z)Vw) = /(x,Z), xe(-l,l), 

(4.19) 

[u(— 1, Z) = Z d _ i, «(1, Z) = Z d . 

The solution is u(x, Z) : [—1, 1] x JR 0 ' — > R. 


4.4 TRADITIONAL NUMERICAL METHODS 

Here we briefly review the traditional methods for solving practical systems with 
random inputs. For the purpose of illustration, we use the simple stochastic ODE 
in example 4.4 as an example. 
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The exact solution to (4.15) is 


u(t, co) = y S(^y)e■“ ( ' u), 


(4.20) 


When the distribution function of a and ft is known, i.e., F„p (a, b) = P(a < a, 
ft < b). the statistical moments of the solution can be evaluated. If a(a>) and ft {co) 
are independent, i.e., F a p(a, b) = F a (a)Fp{b), then 


E [u(t, a>)] = E[ft]E [e- at ] . 


For example, if we further assume ft(co) = 1, that is, the initial condition is not 
random, and a(co) ~ Af (0, 1) is a standard Gaussian random variable with zero 
mean and unit variance, then 



and 


a 2 = E [u 2 ] — (E [m]) 2 = e 21 — e ' . 


4.4.1 Monte Carlo Sampling 

Monte Carlo sampling (MCS) is a statistical sampling method that was popularized 
by physicists from Los Alamos National Laboratory in the United States in the 
1940s. The general procedure for (4.15) is as follows. 

1. Generate identically and independently distributed random numbers Z (,) = 

ft^), i = 1 , ,M, according to the distribution of F a p(a, b). Note 

once again that the dependence structure of a and ft is required to be known. 

2. For each i = 1 solve the governing equation (4.15) and obtain 

u {i \t) = u(t , Z (,) ). 

3. Estimate the required solution statistics. For example, the solution mean can 
be estimated as 



(4.21) 


Other solution statistics can be estimated via proper schemes from the solution en- 
semble {h (,) }. It is obvious that steps 1 and 3 are preprocessing and postprocessing 
steps, respectively. Only step 2 requires solution of the original problem, and it 
involves repetitive simulations of the deterministic counterpart of the problem. 

Error estimate of MCS follows immediately from the Central Limit Theorem 
(theorem 2.26). Since {u(t, Z u> )} are independent and identically distributed ran- 
dom variables, the distribution function of u(t) converges, in the limit of M -» 
oo, to a Gaussian distribution Af(E[u\(t), a 2 (t)/M), whose standard deviation is 
AW l/2 rr„(0, where a u is the standard deviation of the exact solution. Hence the 
widely adopted concept that the error convergence rate of MCS is inversely propor- 
tional to the square root of the number of realizations. 

It is obvious that the MCS procedure can be easily extended to a general system 
such as (4.13). The only requirement is that one needs a well-established solver to 
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solve the corresponding deterministic system. Although very powerful and flexible, 
the convergence rate of MCS, 0(M~ I/2 ), is relatively slow. Generally speaking, if a 
one-digit increase in solution accuracy of the statistics is required, one needs to run 
roughly 100 times more simulations and thus increase the computational effort by 
100 times. For large and complex systems where the solution of a single determinis- 
tic realization is time-consuming, this poses a tremendous numerical challenge. On 
the other hand, a remarkable advantage of MCS lies in the fact that the 0(M~ l/2 ) 
convergence rate is independent of the total number of input random variables. That 
is, the convergence rate is independent of the dimension of the random space. This 
turns out to be an extremely useful property that virtually no other methods possess. 

4.4.2 Moment Equation Approach 

The objective of the moment equation approach is to compute the moments of the 
solution directly because in many cases the moments of the solution are what one 
needs. 

Let p(t) = E [k]; then by taking the expectation of both sides of (4.15) we obtain 

— (t) = — E [an], /r( 0) = E[/3], 

at 

This requires a knowledge of E[aw], which is not known. We then attempt to derive 
an equation for this new quantity. This is done by multiplying (4.15) by a and then 
taking the expectation, 

— E [au\{t) = — E[q! 2 w], E[aw](0) = E[a/3]. 

dt 

A new quantity E [a 2 u] appears. If we attempt a similar approach by multiplying 
the original system by or and then taking the expectation, the equation for the new 
quantity is 

— E [a 2 u](t) = — E[o' 3 m], E[a 2 M](0) = E [a 2 /}], 

dt 

which now requires yet another new quantity, E[a 3 «], In general, if we define 
Jlkit) — E [ot k u\ for k = 0, 1 , ... , then we obtain 

7 '-V-kit ) = -Jlk+ 1, Mir(0) = E[c x k p], 

dt 

The system of equations thus cannot be closed, as it keeps introducing new vari- 
ables that are not included in the existing system. This is the well-known closure 
problem. The typical approach to remedy the difficulty is to assume, for a fixed k, 
that the higher-order term is a function of the lower-order ones, that is, 

fik+i = g(P- o, • • • , JZk), 

where the form of the function g is determined on a problem-dependent basis with 
(hopefully) sufficient justification. There exists, to date, no effective general strat- 
egy for solving the closure problem. Also, the error caused by most, if not all, 
closure techniques is not easy to quantify. 
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4.4.3 Perturbation Method 


In the perturbation method, the random part of the solution is treated as a perturba- 
tion. The fundamental assumption is that such a perturbation is small. And this is 
typically achieved by requiring the standard deviation of the solution to be small. 
To illustrate the idea, let us again use the simple ODE example (4.15). For ease of 
exposition, let us further assume that the mean value of a is zero, i.e., E[ot] = 0, 
and that the initial condition is a fixed value. That is, 
du 

— (f, co) = —a(co)u, u( 0, co) — p. 
dt 

The perturbation method can be applied when the variation of a is small, that is, 
e = 0(a(o))) ~ Oq, 1. If this is the case, we seek to expand the solution as a 
power series of a, 

u(t, co) — Uo(t) + a(co)u\(t) + a 2 (co)u 2 (t ) + • • • , (4.22) 


where the coefficients uq, u\, . . . are supposed to be of same order of magnitude. 
After substituting the expression into the governing equation (4.1), we obtain 

duo du\ ~,du 2 2 

— l-a— h a"— r— H = -a(«o + oiU\ + a 2uo H ). 

dt dt dt 

Since 0(a k ) = e k becomes increasingly small as k increases, we match the terms 
in the equation according to the power of a. 

= 0 , 

= “MO, 

= —Ml, 


0(1 ) : — ^ 


0(e): 7 — 


0(e 2 ): 


duo 

dt 

du\ 

dt 

duo 

dt 


0(e k ): 


duk 

dt 


-Mi-1, 


Similar expansion and term matching of the initial condition result in the initial 
conditions for the coefficients, 

«o(0) = P, Mi(0) = 0, k > 1. 

It is then easy to solve the system recursively and obtain 

t k 

u 0 (t) = P, Ui = -pt, ..., u k = P(-\) k — . 

k\ 

The power series then gives us the solution 

00 —t k 

u(t, co) = V p a k (co), (4.23) 

k\ 

k= 0 

which is the infinite power series of the exact solution u exact (t , co) = p exp(— at). 
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Even though it seems that the perturbation solution can recover the exact solution 
in terms of its power series, the requirement for a to be small is still needed. This 
is because in practice one can use only a finite-term series, which can provide a 
good approximation only when a is small. (Note from (4.23) that when a is 0(1) 
or bigger, the remainder of a finite-term series is always dominant.) 

Derivation of the equations for each terms is done by matching the orders. The 
procedure cannot be easily automated, except for very simple problems such as the 
example here. For practical systems, one usually stops the procedure at the second- 
order terms because of the complexity of the derivation. For first- or second-order 
approximations to be satisfactory, an even stronger requirement for smallness is 
needed. However, a distinct feature of perturbation methods is that, once derived, 
the systems of equations are almost always decoupled and can be solved recur- 
sively. 


Chapter Five 


Generalized Polynomial Chaos 


This chapter is devoted to the fundamental aspects of generalized polynomial chaos 
(gPC). The material is largely based on the work described in [120]. However, the 
exposition here is quite different from the original one in [120], for better under- 
standing of the material. It should also be noted that here we focus on the gPC 
expansion based on globally smooth orthogonal polynomials, in an effort to illus- 
trate the basic ideas, and leave other types of gPC expansion (e.g., those based on 
piecewise polynomials) as research topics outside the scope of this book. 


5.1 DEFINITION IN SINGLE RANDOM VARIABLES 

Let Z be a random variable with a distribution function F z (z ) = P(Z < z) and 
finite moments 

E(|Z| 2m ) = J \z\ 2m dF z (z) < oc, meAf, (5.1) 

where AT = No = [0, 1, . . . } or Af = [0, 1, . . . , N} and for a finite nonnegative 
integer N is an index set. The generalized polynomial chaos basis functions are the 
orthogonal polynomial functions satisfying 

E[4> m (Z)<t>„(Z)] = y n 8 mn , m , n e AT, (5.2) 

where 

Y„ = E[ct> 2 (Z)] , neAf, (5.3) 

are the normalization factors. 

If Z is continuous, then its probability density function (PDF) exists such that 
i dF z (z) = p(z)dz and the orthogonality can be written as 

E[<t>,„(Z)tt>„(Z)] = J <J> H , (z)<b„(z)p(z)dz = y„8 mn , m,nej\f. (5.4) 

Similarly, when Z is discrete, the orthogonality can be written as 

E[cJ>„,(Z)ct>„(Z)] = ^2^ m {z i )^ n (z i )p i = y n 8 mn , m,n eJV. (5.5) 

With a slight abuse of notation, hereafter we will use 

E[/(Z)] = J f(z)dF z (z) 

to include both the continuous case and the discrete case. 
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Obviously, {<J> m (z)} are orthogonal polynomials of z e R with the weight func- 
tion p(z), which is the probability function of the random variable Z. This estab- 
lishes a correspondence between the distribution of the random variable Z and the 
type of orthogonal polynomials of its gPC basis. 


Example 5.1 (Hermite polynomial chaos). Let Z ~ AGO, 1) be a standard Gaus- 
sian random variable with zero mean and unit variance. Its PDF is 


p(z) = 



The orthogonality (5.2) then defines the Hermite orthogonal polynomials {H m (Z)} 
as in (3.19). Therefore, we employ the Hermite polynomials as the basis functions, 

Hq(Z) — 1, Hi(Z) = Z, H 2 (Z) = Z 2 - 1, H 3 (Z) = Z 3 -3Z, .... 


This is the classical Wiener-Hermite polynomial chaos basis ([45]). 


Example 5.2 (Legendre polynomial chaos). Let Z ~ U{— 1 , 1 ) be a random vari- 
able uniformly distributed in (— 1 , 1). Its PDF is p (z) = 1 /2 and is a constant. The 
orthogonality (5.2) then defines the Legendre orthogonal polynomials (3.16), with 

Lo(Z) = l, L\(Z) — Z, L 2 {Z)= 3 -Z 2 - 1 -, .... 

Example 5.3 (Jacobi polynomial chaos). Let Z be a random variable of beta dis- 
tribution in (— 1 , 1 ) with PDF 

p(z) oc (1 — z)“(l + z)P , a, f > 0, 

whose precise definition is in (A. 21). The orthogonality (5.2) then defines the 
Jacobi orthogonal polynomials (A. 20) with the parameters a and f, where 

jt P> (Z) = 1, Jl a ’ P) (Z) = l[ a -^ + (a + p + 2)Z] 

The Legendre polynomial chaos becomes a special case of the Jacobi polynomial 
chaos with a = f = 0. 


In table 5.1, some of the well-known correspondences between the probability 
distribution of Z and its gPC basis polynomials are listed. 


5.1.1 Strong Approximation 

The orthogonality (5.2) ensures that the polynomials can be used as basis functions 
to approximate functions in terms of the random variable Z. 

Definition 5.4 (Strong gPC approximation). Let f(Z) be a function of a random 
variable Z whose probability distribution is Fz(z ) = P(Z < z) and support is 
Iz . A generalized polynomial chaos approximation in a strong sense is fN(Z) e 
Pjv(Z), where P y(Z) is the space of polynomials of Z of degree up to N > 0, such 
that || /(Z) — fsi(Z) || — »• 0 as N oo, in a proper norm defined on Iz- 
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Table 5.1 Correspondence between the Type of Generalized Polynomial Chaos and Their 
Underlying Random Variables" 



Distribution of Z 

gPC basis polynomials 

Support 

Continuous 

Gaussian 

Hermite 

(— 00 , 00 ) 


Gamma 

Laguerre 

[0, 00 ) 


Beta 

Jacobi 

[fl, b] 


Uniform 

Legendre 

[a, b] 

Discrete 

Poisson 

Charlier 

{0,1,2,...} 


Binomial 

Krawtchouk 

{0, 1, 


Negative binomial 

Meixner 

{0,1,2,...} 


Hypergeometric 

Hahn 

{0, 1, 

" N > 0 is a finite integer. 


One obvious strong approximation is the orthogonal projection. Let 

L] Fz (Iz) = { / : 7 z ->• R I nf] < 00 } (5.6) 

be the space of all mean-square integrable functions with norm H/H^ = 
(E[/ 2 ])^ 2 . Then, for any function / e L dFz (Iz), we define its jVth-degree gPC 
orthogonal projection as 

W 1 

Psf = T fk®k(Z), fk = -E[/(Z)<D*(Z)]. (5.7) 

U Yk 

The existence and convergence of the projection follow directly from the classical 
approximation theory; i.e., 

U-Pn/WlI, N-+OC, (5.8) 

dF Z 

which is also often referred to as mean-square convergence. Let P# (Z) be the linear 
space of all polynomials of Z of degree up to A ? ; then the following optimality 
holds: 


\\f ~ Pn fh 2 = inf \\f-g\\ L 2 . 

J nJ " L dF z geV N (Z) 6 " L dF z 


(5.9) 


Though the requirement for convergence (L 2 -integrable) is rather mild, the rate of 
convergence will depend on the smoothness of the function / in terms of Z. The 
smoother / is, the faster the convergence. These results follow immediately from 
the classical results reviewed in chapter 3. 

WhenagPC expansion /w(Z) of a function f(Z) converges to /(Z) in a strong 
norm, such as the mean-square norm of (5.8), it implies that /w(Z) converges to 

P 

f (Z) in probability, i.e., fy -> /, which further implies the convergence in distri- 
bution, i.e., /)v -> /, as ,V 00 . (See the discussion of the modes of convergence 
in section 2.6.) 
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Example 5.5 (Lognormal random variable). Let Y = e x , where X ~ A r (/i. cr 2 ) 
is a Gaussian random variable. The distribution of Y is a lognormal distribution 
whose support is on the nonnegative axis and is widely used in practice to model 
random variables not allowed to have negative values. Its probability density func- 
tion is 

1 (lny-,i > 2 

Py(v) = -j=e 2,2 . (5.10) 

yos/zjt 

To obtain the gPC projection of Y, let Z ~ A/"(0, 1) be the standard Gaussian 
random variable. Then X = n + crZ and Y = f(Z) = e l 'e rT/ . The Hermite 
polynomials should be used because of the Gaussian distribution of Z. By following 
(5.7), we obtain 

Yn(Z ) = e^ ( ° 2/2) J2 ^fiit(Z). (5.11) 

k = 0 K ' 


5.1.2 Weak Approximation 

When approximating a function /(Z) with a gPC expansion that converges 
strongly, e.g., in a mean-square sense, it is necessary to have knowledge of /, that 
is, the explicit form of / in terms of Z. In practice, however, sometimes only the 
probability distribution of / is known. In this case, a gPC expansion in terms of Z 
that converges strongly cannot be constructed because of the lack of information 
about the dependence of f on Z. However, the approximation can still be made to 
converge in a weak sense, e.g., in probability. To be precise, the problem can be 
stated as follows. 

Definition 5.6 (Weak gPC approximation). Let Y be a random variable with dis- 
tribution function Fy(y) — P(Y < y) and let Z be a ( standard ) random variable 
in a set of gPC basis functions. A weak gPC approximation is Yn £ IPjv (Z), where 
Pjv(Z) is the linear space of polynomials of Z of degree up to N > 0, such that Yn 
converges to Y in a weak sense, e.g., in probability. 

Obviously, a strong gPC approximation in definition 5.4 implies a weak ap- 
proximation, not vice versa. We first illustrate the weak approximation via a triv- 
ial example and demonstrate that the gPC weak approximation is not unique. Let 
Y ~ AT//, a 2 ) be a random variable with normal distribution. Naturally we choose 
Z e A/"(0, 1), a standard Gaussian random variable, and the corresponding Hermite 
polynomials as the gPC basis. Then a first-order gPC Hermite expansion 

Ij(Z) = p,Ho + o H\{Z) — p, + cr Z (5.12) 

will have precisely the distribution A/"(/r, a 2 ). Therefore, Y\(Z) can approximate 
the distribution of Y exactly. However, if all that is known is the distribution of Y, 
then one cannot reproduce pathwise realizations of Y via Y\(Z). In fact, I) ( Z ) = 
p 1 1() — a If] (Z) has the same A f (p. , cr 2 ) distribution but entirely different pathwise 
realizations from Y\ . 
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When Y is an arbitrary random variable with only its probability distribution 
known, a direct gPC projection in the form of (5.7) is not possible. More specifi- 
cally, if one seeks an A’th-degi'ee gPC expansion in the form of 

N 

Y N = Y J a k®kW, (5.13) 

k = 0 


with 


at = E[7ct>,(Z)]/ n ., 0 < k < N , (5.14) 

where y k — lE[tt>|] are the normalization factors, then the expectation in the coef- 
ficient evaluation is not properly defined and cannot be carried out, as the depen- 
dence between Y and Z is not known. This was first discussed in [120] where a 
strategy to circumvent the difficulty by using the distribution function F Y (y ) was 
proposed. The resulting gPC expansion turns out to be the weak approximation that 
we defined here. 

By definition, Fy : Iy [0, 1], where Iy is the support of Y . Similarly, F/(z) = 
P(Z < z) : I z — » [0, 1]. Since Fy and Fy map Y and Z, respectively, to a uniform 
distribution in [0, 1], we rewrite the expectation in (5.14) in terms of a uniformly 
distributed random variable in (0, 1], Let U = Fy(Y) = F z (Z) ~ U((), 1); then 
Y = FyUU) and Z = F~\U). (The definition of F~ l is (2.7).) Now (5.14) can 
be rewritten as 

a k = -Eu[Fy\U)^ k (F^(U))] = - [ Ff\u)<$> k {Ff\u))du. (5.15) 

Yk Yk Jo 

This is a properly defined finite integral in [0, 1] and can be evaluated via traditional 
methods (e.g., Gauss quadrature). Here we use the subscript U in Ey to make clear 
that the expectation is over the random variable U . 

Alternatively, one can choose to transform the expectation in (5.14) into the ex- 
pectation in terms of Z by utilizing the fact that Y = Fy 1 (Fz(Z)). Then the ex- 
pectation in (5.14) can be rewritten as 

a k = -¥, z [Fy\F z {Z))^ k (Z)} = — [ Fy l (F z (z))® k (z)dF z (z). (5.16) 
Yk Yk Ji 2 

Though (5.15) and (5.16) take different forms, they are mathematically equivalent. 
The weak convergence of Y N is established in the following result. 

Theorem 5.7. Let Y be a random variable with distribution Fy{y) = P(Y < y) 
and¥i(Y 2 ) < oo. Let Z be a random variable with distribution F z (z) — P(Z < z) 
and let E(|Z| 2m ) < oo for all m e J\f such that its generalized polynomial chaos 
basis functions exist with E, z [<S> m (Z)<t>„(Z)] = 8 mn y n , Vm, n e N . Let 

N 

Y N = J2 a k®k(Z), (5.17) 

k = 0 

where 

a k = -E z [Ff l (F z (Z))® k (Z)l 0 < k < N. 

Yk 


(5.18) 
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Then Y v converges to Y in probability; i.e., 

Y n Y, N oo. (5.19) 

Also, Y tn distribution. 


Proof. Let 

Y = G(Z) = Ff\F z (Z)), 

where G = Ff l oF z : I z — » Iy. Obviously, Y has the same probability distribution 
as that of T, i.e., Fy = Fy, and we have Y = Y and E[T 2 ] < oo. This immediately 
implies 


E[7 2 ] 


f y~ d Fy (y) 

J Iy 

f [Ff l (uj)~du 

Jo 

[ (Ff\F z (z))) 2 dF z {z) 

Jlz 


Therefore, Y = G(Z ) e L 2 dF/ (I z ). Since (5.17) and (5.18) is the orthogonal projec- 
tion of Y by the Mli-degree gPC basis, the strong convergence of Y N to Y in mean 

p ~ ~ p 

square implies convergence in probability, i.e., Y^ —■ ► Y as iV — > oo. Since Y — Y, 
the main conclusion follows. Since convergence in probability implies convergence 
in distribution. Tv -> Y . The completes the proof. 


Example 5.8 (Approximating beta distribution by gPC Hermite expansion). 

Let the probability distribution of T be a beta distribution with probability density 
function p(y) oc (1 — y) a ( 1 + yY . In this case, if one chooses the corresponding 
Jacobi polynomials as the gPC basis function, then the first-order gPC expansion 
can satisfy the distribution exactly. However, suppose one chooses to employ the 
gPC Hermite expansion in terms of Z ~ AM), 1); then a weak approximation can 
still be obtained via the procedure discussed here. All is needed is a numerical ap- 
proximation of the integral (5.15) or (5.16). In figure 5.1, the convergence in PDF 
is shown for different orders of the gPC Hermite expansion. Numerical oscillations 
near the corners of the distributions can be clearly seen, resembling Gibbs oscilla- 
tions. Note that the support of T is in [— 1, 1] and is quite different from the support 
of Z (which is R). 


Example 5.9 (Approximating exponential distribution by gPC Hermite expan- 
sion). Now let us assume that the distribution of T is an exponential distribution 
with the PDF p ( y) oc e~ y . Figure 5.2 shows the convergence of PDF by the gPC 
Hermite expansions. Note that the first-order expansion, which results in a Gaussian 
distribution, is entirely different from the target exponential distribution. As the 
order of expansion is increased, the approximation improves. In this case, if one 
chooses the corresponding gPC basis, i.e., the Laguerre polynomials, then the first- 
order expansion can produce the exponential distribution exactly. 
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Figure 5.1 Approximating beta distributions by gPC Hermite expansions: convergence of 
probability density functions with increasing order of expansions. Left: approxi- 
mation of uniform distribution a = fi = 0. Right: approximation of beta distri- 
bution with a = 2, fi = 0. (More details are in [120].) 



Figure 5.2 Approximating an exponential distribution by gPC Flermite expansions: conver- 
gence of probability density function with increasing order of expansions. 


Example 5.10 (Approximating Gaussian distribution by gPC Jacobi expan- 
sion). Let us assume that the distribution of Y is the standard Gaussian A/"(0, 1) 
and use the gPC Jacobi expansion to approximate the distribution. The conver- 
gence in PDF is shown in figure 5.3, where both the Legendre polynomials and the 
Jacobi polynomials with a — /3 = 10 are used. We observe some numerical oscilla- 
tions when using the gPC Legendre expansion. Again, if we use the corresponding 
gPC basis for Gaussian distribution, the Hermite polynomials, then the first-order 
expansion Y\ — H\(Z) = Z will have precisely the desired .5/(0, 1) distribution. 

It is also worth noting that the approximations by gPC Jacobi chaos with a = 
yS = 10 are quite good, even at the first order. This implies that the beta distribution 
with a = p = 10 is very close to the Gaussian distribution AGO. 1). However, 
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Figure 5.3 Approximations of Gaussian distribution by gPC Jacobi expansions: conver- 
gence of probability density functions with increasing order of expansions. Left: 
approximation by gPC Jacobi polynomials with a = f) = 0 (Legendre polyno- 
mials). Right: approximation by gPC Jacobi polynomials with a = fi = 10. 

a distinct feature of the beta distribution is that it is strictly bounded in a close 
interval. This suggests that in practice when one needs a distribution that is close 
to Gaussian but with strict bounds, mostly because of concerns from a physical 
or mathematical point of view, then the beta distribution can be a good candidate. 
More details on this approximation can be found in appendix B. 

From these examples it is clear that when the corresponding gPC polynomials 
for a given distribution function can be constructed, particularly for the well-known 
cases listed in table 5.1, it is best to use these basis polynomials because a proper 
first-order expansion can produce the given distribution exactly. Using other types 
of polynomials can still result in a convergent series at the cost of inducing approx- 
imation errors and more complex gPC representation. 


5.2 DEFINITION IN MULTIPLE RANDOM VARIABLES 

When more than one independent random variables are involved, multivariate gPC 
expansion is required. Let Z = (Z\, . . . , Z c \) be a random vector with mutually 
independent components and distribution F z (zi , . . . , Zd) — P(Z\ < z\.. . . , Zj < 
z L j). For each i = 1 , . . . , d, let F Zl ( z ,-) = P(Z, < z, ) be the marginal distribution of 
Z, , whose support is I Zj ■ Mutual independence among all Z, implies that F z (z) = 
nti Fz,(.Zi) and I z = I z , x • ■ ■ x I Zd . Also, let {<pk(Zi)}% =0 e P N (Zi) be the 
univariate gPC basis functions in Z, of degree up to N. That is, 

E [(/)„, (Z ,)(!>„ (Zi)] = J 4> m ( Z)(p n (z)dF Zi (z) = S mn y m , 0 <m,n< N. (5.20) 

Let i = (i i, . . . , id) e Nq be a multi-index with |i| = i'i + • • • + id- Then, the 
<:/-variate Ath-degree gPC basis functions are the products of the univariate gPC 
polynomials (5.20) of total degree less than or equal to N\ i.e., 

O i (Z) = 0 il (Z 1 )...^(Z d ), 0< |i| < N. 


(5.21) 
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It follows immediately from (5.20) that 

E [^(ZJOjCZ)] = j ‘Ih (z) 4>j {z)d F z (z) = Vi S u , (5.22) 

where yj = E[<t>?] = y il ■ ■ ■ Yi d are the normalization factors and <5y = <5, Ul • • • 8j d j d 
is the d-\ ariate Kronecker delta function. It is obvious that the span of the polyno- 
mials is F d N , the linear space of all polynomials of degree at most N in d variables. 


K(z) = 



p(Z) = CiW) 


whose dimension is 


dimP^ = 


(V) 


(5.23) 


(5.24) 


The space of homogeneous gPC, following Wiener’s notion of homogeneous chaos, 
is the space spanned by the gPC polynomials in (5.21) of degree precisely N; 
that is, 


VUZ) = 


and 


P ■ h 


dim Vt = 


p{Z ) = 

ih=n 


N + d 
N 


(5.25) 


(5.26) 


The <:/- variate gPC projection follows the univariate projection in a direct manner. 
Let L~ dFz ( Iz ) be the space of all mean-square integrable functions of Z with respect 
to the measure dFz', that is, 




L 


E[/ 2 (Z)] = / j(z)d F z (z) < oo 


(5.27) 


Then for / e L 2 dF (Iz), its A'th-degree gPC orthogonal projection is defined as 


p N /= 


(5.28) 


where 


/t = -E[/<&i] = - [ f (z)<P\(z)d F z (z) , 

Ki Vi J 


V|i| < N. 


The classical approximation theory can be readily applied to obtain 


U-Pn/Wl* 


0, 


N ■ 


and 


(5.29) 


(5.30) 


U-Pn/Wl* = inf \\f-g\\ L > • 


(5.31) 
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Table 5.2 An example of graded lexicographic ordering of the multi-index i in d = 4 di- 
mensions 


i 

Multi-index i 

Single index k 

0 

(0 0 0 0) 

1 

1 

(1 0 0 0) 

2 


(0 10 0) 

3 


(0 0 10) 

4 


(0 0 0 1) 

5 

2 

(2 0 0 0) 

6 


(110 0) 

7 


(10 10) 

8 


(10 0 1) 

9 


(0 2 0 0) 

10 


(0 110) 

11 


(0 10 1) 

12 


(0 0 2 0) 

13 


(0 0 11) 

14 


(0 0 0 2) 

15 

3 

(3 0 0 0) 

16 


(2 10 0) 

17 


(2 0 10) 

18 


Up to this point, the mutual independence among the components of Z has not 
been used explicitly, in the sense that all the expositions are made on dF z (z) and 
Iz in a general manner and do not require the properties of I z = //, x ■ ■ • x I Zd and 
dF z (z ) = dF Zl (z i) • ■ ■ dF Zd (z.d), which are direct consequences of independence. 
This implies that the above presentation of gPC is applicable to more general cases. 

Although clear for the formulation, the multi-index notations employed here are 
cumbersome to manipulate in practice. It is therefore preferable to use a single 
index to express the gPC expansion. A popular choice is the graded lexicographic 
order , where i > j if and only if ]i| > |j | and the first nonzero entry in the difference, 
i — j, is positive. Though other choices exist, the graded lexicographic order is the 
most widely adopted one in practice. The multi-index can now be ordered in an 
ascending order following a single index. For example, for a (d = 4)-dimensional 
case, the graded lexicographic order is shown in table 5.2. 

Let us also remark that the polynomial space (5.23) is not the only choice. An- 
other option is to keep the highest polynomial order for up to N in each direction. 
That is. 


r> N (Z) = 


P-b.-r 


P(Z) = c i®i(Z) 


|l|o<JV 


(5.32) 


where |i|o = maxi< y <^ i r This kind of space is friendly to theoretical analysis (e.g., 
[8]), as properties of one dimension can be more easily extended. On the other hand, 
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dim P'(, = N d . And for large dimensions, the number of basis functions grows too 
fast. Therefore, this space is usually not adopted in practical computations. 


5.3 STATISTICS 

When a sufficiently accurate gPC expansion for a given function f(Z) is avail- 
able, one has in fact an analytical representation of / in terms of Z. Therefore, 
practically all statistical information can be retrieved from the gPC expansion in a 
straightforward manner, either analytically or with minimal computational effort. 

Let us use a random process to illustrate the idea. Consider a process f(t, Z), 
ZeI^ and t e T, where T is an index set. For any fixed t e T , let 

In it, Z)= £ e K 

|1|<JV 

be an ATh-degree gPC approximation of fit , Z); i.e., /’ v ~ / in a proper sense 
(e.g., mean square) for any t e T . Then the mean of / can be approximated as 

X kt)* iU) ) dF z (z) = /o (/), 

|i|£A ) 

(5.33) 

following the orthogonality of the gPC basis functions (5.22). The second moments, 
e.g., the covariance function, can be approximated by, for any t\,t 2 e T, 

c f (t\, t 2 )=mntu z ) - nfitmm, z) -n f m] 

**mMtu Z) - MhMfNiti, z) - /ofe))] 

= X ^kh)kti)l (5.34) 

0<|i|<A r 

The variance of / can be obviously approximated by, for any t e T, 

var(/(f,Z))=E[(/(f,Z)- A t / (f)) 2 ] « X [m^(0] • (5-35) 

0<|i|<A r 

Other statistical quantities of / can also be readily approximated by applying 
their definitions directly to the gPC approximation fy. 


H f (t)±E[f(t,Z)]*E[f N (t,Z)] = 


/( 


Chapter Six 


Stochastic Galerkin Method 


In this chapter we discuss the generalized polynomial chaos (gPC) Galerkin method 
for solving stochastic systems. We first introduce the main idea via a general sto- 
chastic partial differential equation (PDE) and then illustrate more detailed proper- 
ties of the method by applying it to several representative problems. 


6.1 GENERAL PROCEDURE 


Again, without loss of generality, we utilize the stochastic PDE system (4.13). For 
a physical domain D C R £ , l = 1, 2, 3, and T > 0, consider 

u t (x,t,Z) = £(m), Dx(0,r]xI J , 

B{u) = 0, 9Dx[0, T]xR d , (6.1) 

w = uq, D x {t — 0} x R rf , 

where again C is the differential operator, B is the boundary condition operator, uq 
is the initial condition, and Z = (Z\, . . . , Zj) e R rf , d > 1, are a set of mutually 
independent random variables characterizing the random inputs to the governing 
equation. For ease of presentation, let us consider a scalar equation where 

u(x, t, Z ) : D x [0, T] x R rf -* R. 


For a system of equations, the gPC expansion will be applied to each component of 
u individually. 

Let {<J>k(Z)} be the gPC basis functions satisfying 


E[ < D i (Z)d>j(Z)] =Vi (6-2) 

and let P^(Z) be the space of all polynomials of Z e R rf of degree up to A. Then 
the gPC projection of the solution is, for any fixed (x, t), 


N 

Un(x, t , Z) = fq(x, f)<t>i(Z), 
!>l=o 


U\(x, t) = 


— E[m(x, t, Z)<J>j(Z)]. 
Yi 


(6.3) 


Though this is the optimal (in the L 2 dFz sense) approximation in P^, , it is not of 
practical use since the projection requires knowledge of the unknown solution. 

The stochastic Galerkin procedure is a straightforward extension of the classical 
Galerkin approach for deterministic equations. That is, we seek a solution in 
such that the residue of (6.1) is orthogonal to the space K- By utilizing the gPC 
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orthogonal basis functions (6.2), we obtain the following procedure: for any x and t, 
we seek v n € P^ in the form of 

N 

Vn(x, t, Z) = v-,(x, f)<J>i(Z), (6.4) 

|i|=0 

such that for all k satisfying |k| < N, 

Vn(x, t, Z)cJ> k (Z)] = E[£(u JV )'I>k]. D x (0, T], 

E[H(u v )O k ] = 0, dD x [0, T], (6.5) 

Ck = w 0 ,k, D x {t = 0}, 

where Mo, k = E[M ( |<t> k ]/i/ k are the gPC projection coefficients for the initial con- 
dition. Upon evaluating the expectations in (6.5), the dependence in Z disappears. 
The result is a system of (usually coupled) deterministic equations. The size of the 
system is dimP^, = ( N ^ d )- 


6.2 ORDINARY DIFFERENTIAL EQUATIONS 


Let us use the ordinary differential equation in example 4.4 to illustrate the main 
steps of the gPC Galerkin method. 

— ( t , Z) = —a(Z)u, u(t — 0, Z) = /6, 

at 

where the initial condition is assumed to be deterministic (for simplicity). We 
also assume that the random rate constant follows a normal distribution; i.e., a ~ 
Nifi, a 2 ). The corresponding gPC basis will be the Hermite polynomials. Since a 
is the only random input, we need only univariate gPC Hermite expansion 
{Hk{Z)}^ = Q, N > 0, where Z ~ AGO. 1) is the standard normal random variable 
with zero mean and unit variance. The constant a can be expressed as a = /i + aZ. 
Or, in a more general form, 

N 

a N (Z) = J2 a iHi(Z), 

;=o 


where 


no = IX, a\ = cr, a / = 0, i > 1. 


Usually u N is an approximation of a. However, in this case it is an exact expression 
as long as N > 1. Similarly, the initial condition has a trivial gPC projection, 

N 

Av = X>A/,(Z), 

i=0 


where 


bo = P, bt = 0, i > 0, 


which is exact for N > 0. 
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Let 


N 

v N (t,Z) = 

i=0 

be the Mth-degree gPC approximation we seek. The gPC Galerkin procedure 
results in 


E 


d Vn 
dt 


H k 


= E[— a^v^H k ], 


Wk — 0, . . . , N. 


Upon substituting in the gPC expression for a n and v^, we obtain 

. N N 

-? = Y' Y' aiVjdjk Wk = 0, . . . , N, 

dt 


Yk 


( 6 . 6 ) 


where 


e ijk = E [H i (Z)H j (Z)H k (Z)], 0 < i, j, k < N, (6.7) 

are constants. Like the normalization factors y k , these constants can be evaluated 
prior to any computations. In fact, for Hermite polynomials these constants can be 
evaluated analytically, 

y k = k\ k> 0, (6.8) 


i\ j\k\ 

( S - m s-ky: 1 - '• k • and 21 = '■ + > + * is even ' <6 ' 9) 

For other types of gPC basis functions, the analytical expressions for the con- 
stants may not exist. In such cases, one can use numerical quadrature rules with a 
sufficient number of points to compute the constants numerically but exactly since 
the integrands are of polynomial form. 

System (6.6) is thus a system of deterministic ordinary differential equations for 
the coefficients {v k (t)} with initial conditions 

v k (0) = b k , 0 < k < N. (6.10) 

The size of the system is N + 1, and the equations are coupled. Classical numerical 
methods, e.g., Runge-Kutta methods, can be applied, and usually the coupling of 
the system does not pose serious numerical challenges. (More details on numerical 
studies of this problem can be found in [120].) 

We can also rewrite the system in a compact form by using vector notation. By 
taking the summation over i in (6.6) and defining 

1 N 

d jk = Uj Cij k , 

we let A = {Aj k )o<j k <tt be a (N + 1) x (N + 1) matrix. Denote v(t) = 
(vo , . . . , Ojv) r ; then (6.6) can be written as 

= A r v, v(0) = b, (6.11) 


where b = (bo , . . . , bu) T ■ 
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6.3 HYPERBOLIC EQUATIONS 


Let us now consider a simple linear wave equation 

x e (—1, 1), t > 0, 


du(x,t,Z ) 3 u(x,t,Z) 

= c(Z)- 


dt 


dx 


( 6 . 12 ) 


where c(Z) is a random transport velocity that is a function of a random variable 
Z e K. For now we will leave the distribution of Z unspecified and study the 
general properties of the resulting gPC Galerkin system. The initial condition is 
given by 


w(x, 0, Z) = uq(x, Z). 


(6.13) 


The boundary conditions are more complicated, as they depend on the sign of 
the random transport velocity c(Z). A well-posed set of boundary conditions is 
given by 

u(l,t,Z) = u R (t,Z), c(Z) > 0, 

«(— 1, L Z) = ui(t, Z), c(Z) < 0. (6.14) 


The interesting issue to understand is how to properly pose the boundary conditions 
for the gPC Galerkin system. 

Again an univariate gPC expansion is sufficient. For ease of analysis, let us use 
the normalized gPC basis functions, 


E[<J>, (Z)<J>y(Z)] = 8jj, 0 < i, j < N. 


Note that the normalization only requires dividing the nonnormalized basis by the 
square root of the normalization constants. It facilitates the theoretical analysis 
only. In practical implementations, one does not need to normalize the basis. With 
the gPC Galerkin method, we seek, for any (x, t), 


Vn(x, t, Z) = ^ t ),(x, /)0, (Z) 


(6.15) 


and conduct the projection 

dvjf(x, t, Z ) 


E 


dt 




= E 


3uw(x, t, Z) 
c(Z) N \ > ,(Z) 
3x 


for each of the first N + 1 gPC basis k = 0, . . . , N . We obtain 


where 


dVk(x, t ) 

dt 


N 

i=0 


dv, (x, t) 
dx 


k = 0, . . . , N, 


a ik = E[c(Z)d>/(Z)<t> ) t(Z)], 0 <i,k<N. 


(6.16) 


(6.17) 


This is now a coupled system of wave equations of size N + 1 , where the coupling 
is through the random wave speed. If we denote by A the (N + 1) x (N + 1) 
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matrix whose entries are {a,/ i }o<,.k<N, then by definition a,k = aki and A = A 7 is 
symmetric. Let \ — (vq, . . . , v^) T be a vector of length N + 1; then system (6.16) 
can be written as 


3v(x, t ) , 3v(x, t) 

= A 

3 1 3x 


(6.18) 


It is now clear that system (6.18) is symmetric hyperbolic. Therefore, a complete 
set of real eigenvalues and eigenfunctions exists. Moreover, we can understand 
the signs of the eigenvalues, which indicate the direction of the wave for the gPC 
Galerkin system (6. 1 8), based on the signs of wave direction in the original system 
(6.12). 


Theorem 6.1. Consider the gPC Galerkin system (6.18) derived from the original 
system (6.12). Then if c(Z) > 0 (respectively, c(Z) < 0) for all Z, then the eigen- 
values of A. are all nonnegative (respectively, nonpositive); if c(Z ) changes sign, 
i. e. , c(Z) > 0 for some Z and c(Z) < 0 for some other Z, then A has both positive 
and negative eigenvalues for sufficien tly high gPC expansion order N. 

The proof can be found in [48]. 

A less trivial issue is how to impose the inflow-outflow boundary conditions for 
the hyperbolic system (6.18), especially when the wave speed changes signs in the 
original system (6.12). Note that the explicit information about the sign of the wave 
speed disappears in the gPC Galerkin system (6.18). Because (6.18) is symmetric 
hyperbolic, we can diagonalize the system and then impose boundary conditions 
based on the sign of the eigenvalues. 

Since A is symmetric, there exists an orthogonal matrix S r = S' 1 such that 
S r AS = A, where A is a diagonal matrix whose entries are the eigenvalues of A; 
i.e., 

A = diag(A 0 , ...,T j+ , ...,k N ). 

Here the positive eigenvalues occupy indices j = 0 the negatives ones 
occupy indices j = /_,..., N , and the rest, if they exist, are zeros. Obviously, 
7+, j- < N. 

Denote q = (< 70 , ■ ■ ■ , <]n) t = S r v, i.e., qj(x, t ) = X!/t=o s kjVk( x , 0. where sjk 
are the entries for S; then we obtain 


3q(x,f) ^3q(x,f) 


(6.19) 


3 1 3x 

The boundary conditions of this diagonal system are determined by the sign of the 
eigenvalues; i.e., we need to specify 

N 

qj( 1. t) = y^ j s kj u k (\, t), 7 = 0,..., j+, 

k=0 

N 

qj(~ 1, t) = y^ SkjUk(-l, t), 
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Figure 6.1 Convergence property of the gPC Galerkin solution to the wave problem (6.12). 


Left: error convergence with respect to the order of gPC expansion at different 
times. Right: evolution of the solution in mean-square norm in time at different 
gPC orders. 


Here the coefficients u/ c at the boundaries are determined by the exact gPC projec- 
tion of the boundary conditions of u, i.e., u R and u L . Subsequently, the boundary 
conditions for the gPC Galerkin system of equations (6.16) are specified as 


v(l,0 = Sq(l,f), v(— 1, t) = Sq(— 1, t). 


For vanishing eigenvalues, if they exist, no boundary conditions are required. 

It can be shown that the solution of the gPC Galerkin system (6. 1 8) converges to 
the exact solution. In fact, the following error bound was established (see 
theorem 2.2 in [48]), 



(6.20) 


where || • || is the standard L 1 norm in the physical domain (—1, 1), C is a con- 
stant independent of N, t is time, and m > 0 is a real constant depending on the 
smoothness of u in terms of Z. 

A notable feature of the error bound is that the error depends on time in a linear 
manner. This implies that for any fixed gPC expansion order N, the error will grow 
linearly in time. This can be seen in figure 6.1, where the gPC Galerkin solution is 
applied to a simpler version of (6.12) with c(Z) = Z, a uniform random variable in 
(— 1, 1), periodic boundary conditions in space in a domain 0 <x< 2n , and initial 
condition u (x . 0, Z) = cos(x). The exact solution is u ex = cos(x — Zt). On the left 
side of figure 6.1, while we again see the exponential error convergence as the gPC 
order is increased, the error is bigger at a larger time t and requires higher expansion 
orders to reach the converging regime. The time dependence becomes more evident 
on the right side of figure 6.1, where the evolution of the mean-square norm of the 
solutions are plotted. We observe that the gPC Galerkin solutions deviate from the 
exact solution after a certain time. The time at which the accuracy is lost, i.e., errors 
become 0(1), is roughly proportional to the order of the gPC expansion. 
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It is important to realize the following facts. 

• The linear growth of error is not a result of the boundary condition treatment. 
This is obvious from the results in figure 6.1 because they are obtained from 
a problem with a periodic condition (and require no special treatment of the 
boundary conditions). 

• Furthermore, the linear growth of error is not a result of the Galerkin ap- 
proach. In fact, a direct gPC orthogonal projection of the exact solution u ex = 
cos(x — Zt) would require more and more basis function in order to keep 
a fixed projection error. This is because as far as the expansion in Z is con- 
cerned, the time t behaves as a wave number. A larger time t thus requires finer 
representations. This is the fundamental property of approximation theory. 


6.4 DIFFUSION EQUATIONS 

Let us consider a time-dependent stochastic diffusion equation 

^(x, t, Z) = V, • (*(*, Z)V x u) + fix , t), x e D, t e (0, T ], 
at (6.21) 

u{x, 0, Z) = uq(x), u\ao = 0. 

Here we use V v to explicitly specify that the differentiation is in the physical co- 
ordinates x. We assume that the only source of uncertainty is from the diffusivity 
field k, which causes coupling of the resulting gPC Galerkin system of equations. 
Uncertainty in the source term / and initial condition uq will not cause coupling 
and can be dealt with easily. We assume the diffusivity field takes a form 

d 

k(x, Z) = kq(x) + k,(x)Z, , (6.22) 

i = 1 

where if, (x) are deterministic functions obtained by applying a parameterization 
procedure (e.g., the Karhunen-Loeve expansion) to the diffusivity field and that 
Z = (Z \, . . . , Z d ) are mutually independent random variables with specified prob- 
ability distributions. Alternatively, we can write (6.22) as 

d 

k(x, Z) = ki(x)Zj , (6.23) 

i=0 

where Zq = 1 is fixed. For the problem to be well posed, we require 

K(x, Z) > > 0, Vx, Z. (6.24) 

Such a requirement obviously excludes probability distributions of Z that take neg- 
ative values with nonzero probability, e.g., Gaussian distributions. 

Again we seek an Ath-degree gPC approximation 

N 

Vi wit, X, Z) = ^2 £k(f, x)O k (Z), 

|k|=0 
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where E[<J>j(Z)<t>j(Z)] = (5y. Here we again choose to normalize the gPC basis 
first. For ease of exposition, let us adopt the single-index notation discussed in 
section 5.2 and rewrite the gPC expansion as 

A A / N + d\ 

v N (t,x,Z) = J^v i (t,x)0 i (Z), n J, (6.25) 

where the index i is determined by a proper ordering, e.g., the graded lexicographic 
ordering, of the multi-index i e NfJ . 

Upon substituting (6.23) and the gPC expansion (6.25) into the governing equa- 
tion (6.21) and projecting the resulting equation onto the subspace spanned by the 
first M gPC basis polynomials, we obtain, for all A: = 1, . . . , M, 

„ a d M 

• (, Ki(x)S7 x Vj)eij k + f k (t, x ) 

i= 0 7=1 
M 

= V* ' ( a jkWV x Vj) + fk(t,x), (6.26) 

7=1 

where 

e ijk = E[Z,- = J Zi® j(z)® k (z)dF z (z), 0 <i <d, 1 <j,k<M, 


d 

a jk (x) = ’y^Jcj(x)ejjk, 1 < j,k < M, (6.27) 

i=0 


and fk(t, x) are the gPC projection coefficients for the source term /(x, t). (For 
the simple case of deterministic / considered here, f\ — f and fk — 0 for k > 1 .) 

The gPC Galerkin system (6.26) is a coupled system of diffusion equations. It 
can be put into a compact form by using vector matrix notation. Let us denote 
v = (Ci, ... , v m ) t , f = (/i, . . . , /m) 7 . and A(x) = ( a jk )\<j,k<M ■ By definition, 
A = A r is symmetric. The system (6.26) can be written as 


— ( t , x) = V T • [A(x)V x v] + f, (t, x) e (0, T ] x D , 
3 1 

v(0, x) = v 0 (x), v| 3z) = 0, 


(6.28) 


where Vo(x) = (no.i- ■ ■ . , wo ,A/) r is the gPC projection coefficient vector of the 
initial condition «o(x) in (6.21). For the deterministic initial condition considered 
here, mo.i = uq(x) and UQ k = 0 for k > 1 . 

The coupling of the diffusion terms in (6.28) does not pose a serious problem if 
the system is solved explicitly in time. However, an explicit time integration usually 
imposes a severe restriction on the size of the time step because of concerns about 
numerical stability. To circumvent this difficulty, one can employ a semi-implicit 
scheme where the diagonal terms of A are treated implicitly and the off-diagonal 
terms of A are treated explicitly. This results in a naturally uncoupled system to 
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solve with no loss of accuracy in time integration. For example, a first-oder Euler 
forward-backward semi-implicit scheme takes the form 

*7 n + l 


At 


- V x • [D(x)V x V !+1 ] = V x • [S(x)V x V] + f 


(6.29) 


where the superscript n denotes numerical solutions at time level t„, At is the time 
step, and 


D = diag(A), A = D + S. 

Similarly, if we consider the steady-state counterpart of (6.21), 


— V x ■ ( k{x , Z)V A n(x, Z)) = /(x), ieD; 

we find that the gPC Galerkin system is 

u(x, Z)\ dD = 0, 

(6.30) 

-V x ■ [A(x)V x v] = f, x e D; 

v \bd = o. 

(6.31) 


This is a coupled system of elliptic equations. By using the separation of diago- 
nal and off-diagonal terms of A, an efficient iterative scheme can be designed to 
solve the system as an uncoupled system of equations. These algorithms were first 
proposed in [119, 122] and later analyzed in [128]. 


6.5 NONLINEAR PROBLEMS 

The above examples all involve linear problems. This does not imply that the gPC 
Galerkin method can be applied only to linear problems. (In fact, as far as the ran- 
dom space is concerned, none of the examples are linear because the randomness 
in the equations is all in a multiplicative manner.) 

Let us consider the Burgers’ equation from the supersensitivity example in 
section 1.1.1 to illustrate application of the gPC Galerkin method to nonlinear 
problems. 

+ uu x = vu xx , x e r — 1 , 11, 

(6.32) 

[m(-1) = 1 + <$(Z), m(1) = — 1, 

where 8(Z) > 0 is a random perturbation to the left boundary condition at (x = 
— 1) and v > 0 is the viscosity. Again this requires a one-dimensional gPC expan- 
sion. We seek 


Vn(x, t, Z) = ^ Vi(x, t)®i(Z) 


such that 


3 vn 

< 

3 1 


9 V N , 
VN ~; — 
dx 


= vE 


3 2 vn t 

dx 2 


k — 0, .... N. 


By substituting u,y into the equation and using the orthogonality relation of the 
basis functions, we obtain 


dik i v — "\ ^ „ u v j 


ddj 3 2 Vic 

eijk = v - 


dx 2 


k = 0, 


, N, 


(6.33) 
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where e,^ = ]E[ <J>, <I> 7 O*] are constants and yk = E[3>^] are the normalization 
constants (which will be 1 if the basis functions are normalized). 

This is a coupled system of equations where each equation resembles the origi- 
nal Burgers’ equation and the coupling is through the nonlinear term. The classical 
semi-implicit scheme can be applied to solve the system in time, where the non- 
linear coupling terms are treated explicitly and the diffusion terms implicitly. For 
more details, see [123]. 

The nonlinear term uu x in the Burgers’ equation is in quadratic form and results 
in a gPC projection 

= EE*r^*'W 

1=0 7=0 X 

that can be easily evaluated as long as the term is treated as being explicit in time. 
In many cases, however, nonlinear terms in a system do not take polynomial form 
and a direct gPC projection is not straightforward. For example, let us consider the 
projection of a nonlinear term e u , where u is the unknown solution. A gPC Galerkin 
projection requires us to evaluate 

E [e t ' A: O/.] = J e^< i,<s>,(?) <b k {z)dF z {z), (6.34) 

where is the A' th -degree gPC approximation of u. It is clear that the integral over 
z cannot be separated from the summation over i , as in the case of polynomial-type 
nonlinearity. 

A feasible treatment for such kinds of nonlinearity is to approximate the integral 
(6.34) numerically. To this end, one can employ a quadrature rule, or a cubature 
rule in multivariate cases, with sufficient accuracy. That is, 

Q 

E [e VN <$> k ] « J2 (6.35) 

7=1 

where z <JI and w (,} are the nodes and weights of the integration rule in the domain 
defined by the integral. Note that since i>jv(Z) takes a known polynomial form, the 
evaluation of e VN at any node is a simple exercise in polynomial evaluation. 


dvj v , 

vjf — — 
ox 


Chapter Seven 


Stochastic Collocation Method 


In this chapter we discuss the basic ideas behind the stochastic collocation (SC) 
method, also referred to as the probabilistic collocation method (PCM). The collo- 
cation methods are a popular choice for complex systems where well-established 
deterministic codes exist. We first clarify the notion of stochastic collocation, for 
the purposes of this book, and then discuss the major numerical approaches. As in 
the rest of the book, we discuss only the fundamental aspects of SC and leave the 
more advanced research issues untouched. This is particularly true for this chapter 
because SC has undergone rapid development after its systematic introduction in 
[ 118 ]. 


7.1 DEFINITION AND GENERAL PROCEDURE 

In deterministic numerical analysis, collocation methods are those that require the 
residue of the governing equations to be zero at discreet nodes in the computa- 
tional domain. The nodes are called collocation points. The same definition can 
be extended to stochastic simulations. Let us use the stochastic partial differential 
equation (PDE) system (4.13) again to explain the idea, 

u t {x, t, Z) = C(u), D x (0, T ] x Iz , 

B(u) = 0, 3D x [0, T] x I z , (7.1) 

u = mo, D x {/ = 0} x Iz, 

where I 7 C VV , d > 1, is the support of Z. For any given x and t, let )/;(•, Z) be 
a numerical approximation of u. In general, w(-, Z ) ~ m(-, Z) in a proper sense in 
Iz, and the system (7.1) cannot be satisfied for all Z after substituting u for w. 

Let ©m = {Z (/, }^ =| c Iz be a set of (prescribed) nodes in the random space, 
where M > 1 is the number of nodes. Then in the collocation method, for all 
j = 1, . . . , M, we enforce (7.1) at the node Z (/i by solving 

u t (x, t, ZCl) = C{u), D x (0, T], 

B{u) = 0, 3D x [0, T], (7.2) 

u = Mo, D x [t = 0). 

It is easy to see that for each j, (7.2) is a deterministic problem because the value of 
the random parameter Z is fixed. Therefore, solving the system poses no difficulty 
provided one has a well-established deterministic algorithm. Let m (/) = n(-, Z {jy ), 
j = I .... , M , he the solution of the above problem. The result of solving (7.2) is 
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an ensemble of deterministic solutions [u^}J =l . And one can apply various post- 
processing operations to the ensemble to extract useful information about u(Z). 

From this point of view, all classical sampling methods belong to the class of 
collocation methods. For example, in Monte Carlo sampling, the nodal set ©m is 
generated randomly according to the distribution of Z, and the ensemble average 
is used to estimate the solution statistics, e.g., mean and variance. In determinis- 
tic sampling methods, the nodal set is typically the nodes of a cubature rule (i.e., 
quadrature rule in multidimensional space) defined on I z such that one can use 
the integration rule defined by the cubature to estimate the solution statistics. Con- 
vergence of these classical sampling methods is then based on the convergence of 
solution statistics, e.g., moments, resulting in convergence in a weak measure such 
as convergence in distribution. 

In this book we do not label the classical sampling methods as stochastic collo- 
cation. Instead we reserve the term “stochastic collocation” for the type of collo- 
cation methods that result in a strong convergence, e.g., mean-square convergence, 
to the true solution. This is typically achieved by utilizing the classical multivari- 
ate approximation theory to strategically locate the collocation nodes to construct 
a polynomial approximation to the solution. 

Definition 7.1 (Stochastic collocation). Let @ M — {Z (/) | / v ^ | c Iz be a set of 

(prescribed) nodes in the random space, where M > I is the number of nodes, 
and let be the solution of the governing equation (7.2). Then find w(Z ) e 

n(Z) in a proper polynomial space F1(Z) such that w(Z) is an approximation to 
the true solution u(Z) in the sense that ||ic(Z) — u(Z ) || is sufficiently small in a 
strong norm defined on Iz- 

Convergence of stochastic collocation thus requires 

||to(Z) - u(Z)\\ -* 0, M-»o o, 

where the norm is to be determined and is typically an L p norm. 

As of the writing of this book, there exist two major approaches for high-order 
stochastic collocation: the interpolation approach and the discrete projection ap- 
proach (the pseudospectral approach). 


7.2 INTERPOLATION APPROACH 

Interpolation is a natural approach to the stochastic collocation problem defined 
in definition 7.1. The problem can now be posed as follows: given the nodal set 
&m C Iz and {n^)}yLj, find a polynomial w(Z) e Tl(Z) such that yu(Z { ^) = id J> 
for all 7 = 1,..., M. 

The goal can be easily accomplished, at least in principle. One way is to use a 
Lagrange interpolation approach. That is, let 

M 

w(Z ) = J2 l ‘(Z U> )Lj(Z), 
j= i 


(7.3) 
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where 

Lj(Z {,) ) = Stj, 1 < i, j < M, (7.4) 

are the Lagrange interpolating polynomials. The approach, albeit straightforward 
in formulation, can become nontrivial in practice. This is mostly due to the fact 
that unlike the situation in univariate interpolation, where ample mathematical the- 
ory exists, many fundamental issues of multivariate Lagrange interpolation (when 
d > 1) are not clear. Issues such as the existence of Lagrange interpolating polyno- 
mials for any given set of notes are not well understood. 

The other way is a matrix inversion approach , where we prescribe the polyno- 
mial interpolating basis first. For example, let us use a set of gPC polynomial bases 
<I>k(Z) and construct 

N 

w N (Z) = ^2 u>k<F k (Z) 

|k|=0 

as the gPC approximation of u(Z). The interpolation condition w(Z iJ] ) = u (j> 
results in the following problem for the unknown expansion coefficients 

A r w = u, 

where 

A = (ch k (Z <7) ), 0 < |k| < N, 1 < j < M, 

is the Vandermonde-like coefficient matrix, w is the vector of the expansion coef- 
ficients, and u = (m(Z (1) ), . . . , u(Z lM> )) T . To prevent the problem from becoming 
underdetermined, we require the number of collocation points not to be smaller 
than the number of gPC expansion terms, i.e., M > ( N t d )- The advantage of the 
matrix inversion approach is that the interpolating polynomials are prescribed and 
well defined. Once the nodal set is given, the existence of the interpolation can 
always be determined in the spirit of determining whether the determinant of A 
is zero. However, an important and very practical concern is the accuracy of the 
interpolation. Even though the interpolation has no error at the nodal points, error 
can become wild between the nodes. This is particularly true in high-dimensional 
spaces. Here again, we find rigorous analysis lacking and no general (and sound) 
guideline for choosing the location of the nodes. Many ad hoc choices do exist, for 
example, those based on design of experiments (DoE) principles. However, none 
has become satisfactory for general purposes. 

Since univariate interpolation is a well-studied topic, one solution to multivari- 
ate interpolation is to employ a univariate interpolation and then fill up the entire 
space dimension by dimension. By doing so the properties and error estimates of 
univariate interpolation can be retained as much as possible. In fact, the afore- 
mentioned two approaches, the Lagrange interpolation and matrix inversion ap- 
proaches, are direct conceptual extensions of the univariate interpolation techniques 
in section 3.4. Let us recall that in the univariate case d = 1, i.e., Z e R. Let 
(Z (1 \ . . . , Z^+T) be a set of distinct nodes and let (n (1) , . . . , u ( - N+1> ) be the so- 
lution at the nodes. Then an interpolation polynomial n n/(Z) that interpolates 
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a given function /(Z) can be constructed either in the Lagrange form (3.42) or 
by inverting the Vandermonde matrix to obtain (3.43). Note that the two are equi- 
valent because of the uniqueness of univariate interpolation. It is also understood 
that the interpolating nodes offering high accuracy are the zeros of the orthogonal 
polynomials {<t>,t(Z)}. 

7.2.1 Tensor Product Collocation 

For multivariate cases with d > 1, for any 1 < i < d, let Q mj be an interpolating 
operator such that 

QmXf] = n e p m ,(z,) 

is an interpolating polynomial of degree m , , for a given function / in the Z, variable 
by using m t + 1 distinct nodes in the set ©'”' = {Z ( . , . . . , Z ( <m,) }. Then the most 
straightforward approach to interpolating / in the entire space Iz C is to use a 
tensor product approach. That is, 

8m = Qm, ® ® Qm d , (7.5) 

and the nodal set is 

© M = x x ©“■*, (7.6) 

where the total number of nodes is M = m\ x x mj. 

By using the tensor product construction, all the properties of the underlying 
one-dimensional interpolation scheme can be retained. And error estimate in the 
entire space can be easily derived. For example, let us assume that the number of 
points in each dimension is a constant; i.e., m \ = ■ ■ ■ = m d = m, and that the 
one-dimensional interpolation error in each dimension 1 < i < d follows 

(I- &»,)[/] cx m-“, 

where the constant a > 0 depends on the smoothness of the function /. Then the 
overall interpolation error also follows the same convergence rate 

(/ - 2 m )[./1 oc m~ a . 

However, if we measure the convergence in terms of the total number of points, 
M — m d in this case, then 

(/-0 M )[/]aM-^, d> 1. 

For large dimensions d 1 , the rate of convergence deteriorates drastically and 
we observe very slow convergence, if there is any, in terms of the total number of 
collocation points. Moreover, the total number of points, 

M = m d , 

grows very fast for large d. This poses a numerical challenge because each collo- 
cation point requires a simulation of the full-scale underlying deterministic system, 
which can be time-consuming. This is the well-known curse of dimensionality. 
For this reason, tensor product construction is mostly used for low-dimensional 
problems with d typically less than 5. A detailed theoretical analysis for stochastic 
diffusion equations can be found in [7]. 
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7.2.2 Sparse Grid Collocation 

An alternative approach is Smolyak sparse grids. A detailed presentation of the 
Smolyak construction, originally proposed in [96], is beyond the scope of this 
entry-level textbook, and we refer interested readers to the many more recent stud- 
ies. It is sufficient, for the purposes of this book, to know that the Smolyak sparse 
grids are still based on tensor product construction but are only a subset of the full 
tensor grids. The construction takes the following form ([114]): 

e*= E (7-7) 

N- </+l<|i|<JV ' 1 

where N > d is an integer denoting the level of the construction. Though the 
expression is rather complex, (7.7) is nevertheless a combination of the subsets of 
the full tensor grids. The nodal set, the sparse grids, is 

©m = U (0? x ■ ■ ■ x 0'“). (7.8) 

A"— <7+l<|i|<A/ 

Again it is clear that this is the union of a collection of subsets of the full ten- 
sor grids. Unfortunately, there is usually no explicit formula to determine the total 
number of nodes M in terms of d and N. 

Because the construction in (7.7) employs one-dimensional interpolations with 
various numbers of nodes, it is preferable that one-dimensional nodal sets be nested. 
That is, the one-dimensional nodal sets satisfy 

©i C ©{, i < j. (7.9) 

If this condition is met, then the total number of nodes in (7.8) can reach a mini- 
mum. However, in practice, since one-dimensional nodes are typically the zeros of 
orthogonal polynomials, the nested condition (7.9) is usually not satisfied. 

One popular choice of nested grids is Clenshaw-Curtis nodes, which are the ex- 
trema of Chebyshev polynomials and are defined as, for any 1 < i < d, 

z\ j) = - cos ~ l) , j=l,...,m-, (7.10) 

ntj — 1 

where an additional index is introduced via the superscript k. With a slight abuse 
of notation, we will use k to index the point sets instead of using the total num- 
ber of points w,. Let the point sets double with an increasing index of k > 1, i.e., 
= 2 k ~ 1 + 1, and define m\ = 1 and Z ( (l) = 0. It is easy to see that because 
of the doubling of the nodes we have ©j C ©i +1 and that the sets are nested. 
The additional index k here is often referred to as the level of Clenshaw-Curtis 
grids. The higher the level, the finer the grids. For a more detailed discussion of 
Clenshaw-Curtis nodes, see [27], 

By using Clenshaw-Curtis grids as one-dimensional nodes, the Smolyak con- 
struction (7.7) can be expressed in terms of the level k as well. Let N = d + k, 
where k > 0, and then the “nestedness” of the base one-dimensional nodes can 
be retained, 0 k c ©i+i- Again here we do not use the total number of nodes M, 
which does not have an explicit and closed form expression in terms of d and k, 
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Figure 7.1 Two-dimensional ( d = 2) nodes based on the same one-dimensional extrema of 
Chebyshev polynomials at level k = 5. Left: tensor grids. The total number of 
points is 1089. Right: Smolyak sparse grids. The total number of nodes is 145. 


to index the sets. In the more interesting case of high-dimensional spaces, the total 
number of points satisfies the following estimate: 

M = #@a- ~ 2 k d k /k\, c/ » 1. (7.11) 

It has been shown ([9]) that interpolation through the Clenshaw-Curtis sparse grid 
interpolation is exact if the function is in P^. (In fact, the polynomial space for 
which the interpolation is exact is slightly bigger than Pj?.) For large dimensions 
d y>> 1, dim IP^ = ( d ^ k ) ~ d k /k\. Therefore, the number of points from (7.11) is 
about 2 k more and the factor is independent of the dimension d. For this reason, the 
Clenshaw-Curtis-based sparse grids construction is sometimes regarded as optimal 
in high dimensions. There have been extensive studies on the approximation prop- 
erties of sparse grids, particularly those based on Clenshaw-Curtis nodes. Here we 
cite one of the early results from [9]. For functions in space F d = {f : [— 1, 1 ] £/ -* 
Rl^' 1 / continuous, ij < l, V j), the interpolation error follows 

WI-QmU < C d , e M~ e (log 

where M is the total number of nodes. Compared to that for tensor grids, the curse 
of dimensionality, albeit still present, is lessened. 

An example of two-dimensional sparse grids is shown in figure 7.1, where we 
observe a significant reduction in the number of nodes. 


7.3 DISCRETE PROJECTION: PSEUDOSPECTRAL APPROACH 

Another approach to achieving the goal of stochastic collocation, as defined in 
definition 7.1, is to conduct discrete projection, or the pseudospectral approach 
(as it is termed in [116]). To this end, let us first recall the notion of the quadrature 
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rule in section 3.5 and extend it to multivariate space. Often termed the cubature 
rule, it is an integration rule that seeks to approximate an integral 

J f(z)dF z (z), zeM rf , d> 1, (7.12) 

by 

Q 

UC[f\=J2f( zU) ) a(i) ' e>!, (7.13) 

7=1 

where (z^\ a ( 7*), j = I Q, are the nodes and their corresponding weights. 

The integration rule is convergent if it converges to the integral (7.12) as Q -> oo. 
Typically, the accuracy of an integration is measured by polynomial exactness. An 
integration rule of degree m implies that the approximation (7.13) is exact for any 
integrand / that is a polynomial of degree up to m and is not exact for at least one 
polynomial of degree m + 1 . Hereafter we will freely interchange the notation of 
the integration rule and the cubature rule, with an understanding that in univariate 
cases it is reduced to the quadrature rule. 

To conduct discrete gPC projection, we recall the exact orthogonal gPC projec- 
tion of m(Z), the solution of (7.1), 

N 

u N (Z) = P N u = « k d> k (Z), (7.14) 

|k|=0 

where the expansion coefficients are obtained as 

«k = -E[«(Z)O k (Z)] = — [ u(z)<b k (z)dF z (z), V|k| < N, (7.15) 
Kk Yk J 

where y k = E[<t>£] are the normalization constants of the basis. 

The idea of discrete projection is to approximate the integrals in the expan- 
sion coefficients (7.15) of the continuous generalized polynomial chaos (gPC) pro- 
jection (7.14) by an integration rule. The discrete projection of the solution of 
(7.1) is 

N 

w n (Z)=J2 {u ^(Z), (7.16) 

|k|=0 

where the expansion coefficients are 

1 1 ^ 

w k = — U Q [u(Z)<S> k (Z)} = -y«(z w )$ k (z w )a w . (7.17) 

Yk Ykyr[ 

It is clear that by using the cubature rule U - the coefficients {u) k } are approxima- 
tions to the exact projection coefficients {z) k } in (7.15). Subsequently, the discrete 
projection wn(Z) approximates the continuous projection u ,\ ( Z ) of (7.14). More- 
over, if the cubature rule is convergent, then w k converges to u k as Q —> oo, and 
u>n and u ,y become identical. The following result is then straightforward. 
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Proposition 7.2. For u(Z) e L dFz (Iz), let u ,\ ( Z ) be the gPC projection defined in 
(7.14) and (7.15) and let w^(Z) be the discrete gPC projection defined in (7.16) 
and (7.17). Assume that the cubature rule U~ used in (7.17) is convergent; then as 
Q oo, u>k — »■ iik, for all |k| < N, and 


wpj(Z) -> un(Z), VZ. 


(7.18) 


The error induced by can be easily separated by the triangular inequality. 

Proposition 7.3. For u(Z) e L dFz (Iz), let u ,\ ( Z ) be the gPC projection defined in 
(7.14) and (7.15) and let w^(Z) be the discrete gPC projection defined in (7.16) 
and (7.17). Then, 



The first term of the error is the gPC projection error induced by using finite-order 
(ATh-order) polynomials. The second term of the error is the difference between 
the continuous gPC projection and the discrete projection and is caused by using a 
cubature rule with finite accuracy. It can be expressed as 



(7.20) 


and is termed “aliasing error” in [ 1 16], by following the similar nomenclature from 
classical deterministic spectral methods (cf. [13, 46]). 

Similar to the interpolation approach, the construction of gPC expansion via 
(7.16) and (7.17) can also be considered a postprocessing step after all the com- 
putations are finished at the cubature nodes. A distinct feature of the discrete pro- 
jection approach is that one can compute only the coefficients that are important 
for a given problem without evaluating the rest of the expansion coefficients. This 
may happen, for example, when global sensitivity is required for some input ran- 
dom variables Z. This is in contrast to the gPC Galerkin method, where all the gPC 
coefficients are coupled and solved simultaneously. 

Since the main task in the discrete gPC projection approach is to approximate the 
integrals in (7.15), the problem of multivariate polynomial approximation is trans- 
formed to a problem of multivariate integration, where the accuracy of the chosen 
integration rule is critical. Compared to multivariate interpolation, which is used 
by the stochastic interpolation collocation approach, there exist, relatively speak- 
ing, more results on multivariate integration, which is nevertheless a challenging 
and very active research topic. 

7.3.1 Structured Nodes: Tensor and Sparse Tensor Constructions 

Since Gauss quadrature rules offer high accuracy for univariate integrations, it is 
natural to extend them to multivariate integrations. The most straightforward way of 
constructing high-order integration rules is to extend quadrature rules (in univariate 
cases) to high-dimensional spaces by using tensor construction. Let U m ‘ be a Gauss 
quadrature rule in the Z, direction of Z = (Z i, . . . , Z £ /), d > 1, with a nodal 
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set 0'"' consisting of an m, > 1 number of nodes. Let us assume that it is exact for 
all polynomials in P 2 ” 1 ' -1 (Z,). Then a tensor construction is 

U Q = U mi 

Obviously, the nodal set is 

® Q = x ••• x @'" d , 

whose total number of nodes is Q = m\ x ■ ■ • x /«,/. This integration rule is exact 
for all polynomials in 

P 2 ""— 1 (Zi)<g>---<g>P 2m ‘'- 1 (Z rf ). 

Though easy to construct and of high accuracy, the problem is again the rapid 
growth of the total number of points in high-dimensional random spaces. If we 
use an equal number of nodes in all directions, nt\ = ■ ■ ■ = m.d = m, then the total 
number of nodes is Q — m d . For d 1, this can be a staggeringly large number. 
(Again let us keep in mind that at each node the full-scale deterministic problem 
needs to be solved.) Consequently, the tensor product approach is mostly used at 
lower dimensions, e.g., d < 5. 

To reduce the total number of nodes while keeping most of the high accuracy of- 
fered by Gauss quadrature, the Smolyak sparse grids construction can be employed, 
similarly to the sparse interpolation discussed in section 7.2.2, 

U Q = J2 (-l) J ' Hm| -( V (U mi ®---®U md ), (7.21) 

where N > d is an integer denoting the level of construction. The grid set, the 
sparse grids, is 

© e = (J (©7' x • • ■ x ©70- (7.22) 

jV-rf+l<|m|<V 

Again it is clear that this is the union of a collection of subsets of the full ten- 
sor grids. Usually there is no closed-form explicit formula for the total number of 
nodes Q. 

Depending on the choice of Gauss quadrature in one dimension, there are a 
variety of sparse grid constructions. And they offer different accuracy. Many of 
the constructions are based on the Clenshaw-Curtis rule in one dimension, whose 
properties are closely examined in [108]. Here we will not engage in further in- 
depth discussion of the technical details. Interested readers should see, for example, 
[9, 39, 83, 84], 

7.3.2 Nonstructured Nodes: Cubature 

The study of cubature rules is a relatively old topic but is still actively pursued. The 
goal is to construct a cubature rule with a high degree of polynomial exactness and a 
lesser number of nodes. To achieve this, structured nodes are usually not considered 
and many studies are based on geometric considerations. Most of the rules are given 
in explicit formulas, regarding the node locations and their corresponding weights. 
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Depending on the required accuracy for the discrete gPC projection, one can choose 
a proper cubature with an affordable number of nodes for simulations. For extensive 
reviews and collections of available cubature rules, see, for example, [21, 49, 99]. 

It is worth pointing out that the classical error estimate, in terms of error bounds 
and such, is nontrivial to carry out for cubature rules. Consequently, the accuracy 
of cubature rules is almost always classified by their polynomial exactness. 


7.4 DISCUSSION: GALERKIN VERSUS COLLOCATION 

While the gPC expansion provides a solid framework for stochastic computations, 
a natural question to ask is. For a given practical stochastic system, should one use 
the Galerkin method or the collocation method? 

The advantage of stochastic collocation is clear — ease of implementation. The 
algorithms are straightforward: (1) choose a set of nodes according to either multi- 
variate interpolation theory or integration theory; (2) run deterministic code at each 
node; and (3) postprocess to construct the gPC polynomials, via either the interpo- 
lation approach (section 7.2) or the discrete projection approach (section 7.3). The 
applicability of stochastic collocation is not affected by the complexity or nonlin- 
earity of the original problem so long as one can develop a reliable deterministic 
solver. The executions of the deterministic algorithm at each node are completely 
independent of each other and embarrassingly parallel. For these reasons, stochastic 
collocation methods have become very popular. 

The stochastic Galerkin method, on the other hand, is more cumbersome to im- 
plement. The Galerkin system (6.5) needs to be derived, and the resulting equations 
for the expansion coefficients are almost always coupled. Hence new codes need to 
be developed to deal with the larger, coupled system of equations. When the origi- 
nal problem (6.1) takes highly complex and nonlinear form, the explicit derivation 
of the gPC equations can be nontrivial — sometimes impossible. 

However, an important issue to keep in mind is the accuracy of the methods. The 
stochastic Galerkin approach ensures that the residue of the stochastic governing 
equations is orthogonal to the linear space spanned by the gPC polynomials, as in 
(6.5). In this sense, the accuracy is optimal. On the other hand, stochastic colloca- 
tion approaches, with no error at the nodes, introduce errors either because of the 
interpolation scheme (if interpolation collocation is used) or because of the inte- 
gration rule (if discrete projection collocation is used). Both errors are caused by 
introduction of the nodal sets and can be classified as aliasing errors. Though in one 
dimension, aliasing error can be kept at the same order as the error of the finite order 
Galerkin method, in multidimensional spaces the aliasing errors can be much more 
significant. Roughly speaking, at a fixed accuracy, which is usually measured in 
terms of the polynomial exactness of the approximation, all of the existing colloca- 
tion methods require the solution of a (much) larger number of equations than that 
required by the gPC Galerkin method, especially for higher-dimensional random 
spaces. This suggests that the gPC Galerkin method offers the most accurate solu- 
tions involving the least number of equations in multidimensional random spaces, 
even though the equations are coupled. 
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The exact cost comparison between the Galerkin and the collocation methods 
depends on many factors including error analysis, which is largely unknown, and 
even on coding efforts involved in developing a stochastic Galerkin code. How- 
ever, it is safe to state that for large-scale simulations where a single determin- 
istic computation is time-consuming, the stochastic Galerkin method should be 
preferred (because of the lesser number of equations) whenever (1) the coupling 
of gPC Galerkin equations does not incur much additional computational and cod- 
ing effort. For example, for Navier-Stokes equations with random boundary/initial 
conditions the evaluations of the coupling terms are trivial ([121]) or (2) efficient 
solvers can be developed to effectively decouple the gPC Galerkin system. For ex- 
ample, for stochastic diffusion equations decoupling can be achieved in the manner 
in (6.29). 

Another factor that should be taken into account as part of the effort in imple- 
menting the stochastic Galerkin method is that the properties of the Galerkin system 
may not be clear, even when the baseline deterministic system is well understood. 
And this may affect our design of numerical algorithms for the Galerkin system. A 
simple example is the linear wave equation with random wave speed in section 6.3. 


Chapter Eight 


Miscellaneous Topics and Applications 


The discussions up to this point have involved numerical algorithms, mostly based 
on generalized polynomial chaos (gPC), for propagating uncertainty from inputs to 
outputs. Here we will discuss several related topics regarding variations and appli- 
cations of gPC methods. More specifically, we will consider efficient algorithms 
for 

• random geometry, where the uncertain input is in the specification of the com- 
putational domain, 

• parameter estimation, i.e.. how to estimate the probability distribution of the 
input parameters, and 

• uncertainty in models and how to “correct” it by using available measurement 
data. 

For the second and third topics, the help of experimental measurement is required. 
gPC -based stochastic methods can improve the accuracy of existing methods at 
virtually no simulation cost. 

Unlike those in the previous chapters of this book, the topics in this chapter are 
closer to research issues. However, here we will only briefly discuss these topics 
with a focus on their direct connections with gPC methods. More in-depth general 
discussions and literature reviews can be found in the references. 


8.1 RANDOM DOMAIN PROBLEM 

Throughout this book, we have always assumed that the computational domain is 
fixed and contains no uncertainty. In practice, however, it can be a major source 
of uncertainty, as in many applications the physical domain cannot be determined 
precisely. The problem with uncertain geometry, e.g., a rough boundary, has been 
studied in areas such as wave scattering with many specially designed techniques. 
(See, for example, a review in [113].) For general-purpose partial differential equa- 
tions (PDEs), however, numerical techniques in uncertain domains are less devel- 
oped. Here we discuss general numerical approaches by following one of the earlier 
systematic studies in [130]. 

Let D(a>), a> e Q, be a random domain whose boundary dD(co) is the source of 
randomness. Since dD(co) is a random process, we first seek to parameterize it by a 
function of a finite number of independent random variables. This parameterization 
procedure is required for the ensuing stochastic simulations. Its implementation, 
though more or less straightforward on a conceptual level, can be nontrivial in 
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practice. (Detailed discussions are in section 4.2.) Let 3 D(Z) ,ZeR d ,d> 1, 
be the parameterization of the random boundary. A partial differential equation 
defined on this domain can be written as 

u t {x , t ) = £(x; u), D(Z) x (0, £], 

B(u) = 0, 3Z)(Z) x [0, T], (8.1) 

u = wo, D(Z) x {t = 0}, 

where x = (xi, . . . , Xf ) , i = 1, 2, 3, is the coordinate in the random domain D(Z). 
For simplicity, here the only source of uncertainty is assumed to be from 3 D(Z). 
Note that even though the form of the governing equations is deterministic (it does 
not need to be), the solution still depends on the random variables Z and is a sto- 
chastic quantity. That is, 

w(x, t ) : D(Z) x [0, T ] — > R”“ 

depends implicitly on the random variables Z e 

The key to solving (8.1) is to define the problem on a fixed domain where the op- 
erations for statistical averaging become meaningful. A general approach is to use 
a one-to-one mapping ([130]). Let y — (yi , . . . , y f ) , l = 1,2, 3, be the coordinate 
in a fixed domain fcR ! and let 

y = y(x,Z), x = x(_y, Z), VZ e IT', (8.2) 

be a one-to-one mapping and its inverse such that the random domain D(Z) can be 
uniquely transformed to the deterministic domain E. Then (8.1) can be transformed 
to the following: for all Z e find v = v(y, Z) : E x -> K"“ such that 

v t (y, t, Z) — L(y, Z; v), Ex( 0,T]xR d , 

B(v ) = 0, 3 Ex [0, T] x M d , (8.3) 

v = vo, E x {t = 0} x 

where the operators L and B are transformed from £ and B, respectively, and ty 
is transformed from uq because of the random mapping (8.2). The transformed 
problem (8.3) is a stochastic PDE in a fixed domain, and the standard numerical 
techniques, including those based on gPC methodology, can be readily applied. 

The mapping technique seeks to transform a problem defined in a random do- 
main into a stochastic problem defined in a fixed domain. The randomness in the 
domain specification is absorbed into the mapping and further into the definition 
of the transformed equation. Thus, it is crucial to construct a unique and invertible 
mapping that is also robust and efficient in practice. For some domains, this can be 
achieved analytically ([102]). 

Example 8.1 (Mapping for a random channel domain). Consider a straight 
channel in two dimensions, with length L and height H. Let us assume that the 
bottom boundary is a random process with zero mean value and other known dis- 
tribution functions. That is, the channel is defined as 


(xi, xo) e D(a>) = [0, L] x [h(x , co), H], 


(8.4) 
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where E[/z(jc, «)] = 0. It is easy to see that a simple mapping 

H 


y i = *i 

can map the domain into 


yi 


H — h(x, co) 
(y u y 2 )eE = [0,L]x[0, H]. 


(x 2 — h(x, co)) 


Example 8.2 (Mapping of a diffusion equation). Consider a deterministic Pois- 
son’s equation with homogeneous Dirichlet boundary conditions in a random do- 
main D(Z(w)), Z e Iz, 


V ■ [c(x)Vh(x, Z)]=a(x) in D{Z), 
u(x,Z)=0 on 3 D(Z), 


(8.5) 


where no randomness exists in the diffusivity field c(x) and the source field a (x ) . 
The stochastic mapping (8.2) transforms (8.5) into a stochastic Poisson’s equation 
in a deterministic domain E : 


E 

i=i 


3 

3 yt 


K 



<Xij(y, Z) 


dv \ 

Byj) 


J~ 1 f(y, Z) in Exl z , (8.6) 


v(y, Z) — 0 on dE x Iz, 


where the random fields k and / are the transformations of c and a, respectively, 
J is the transformation Jacobian 


J(y, Z) = 


3(.vt — ,yi) 
d(xi ,xt)' 


and 


a t j {y, Z) = J 1 V y; • V yj , 1 < i, j < i. (8.7) 

Though (8.6) is more complex, it is a stochastic diffusion problem in a fixed do- 
main. The existing methods, such as those based on gPC, can be readily applied. 


Example 8.3 (Diffusion in a random channel domain). Now let us combine the 
aforementioned two examples and consider diffusion problem in a random channel 
domain. This is the same example as that used in [130]. 

Consider the steady-state diffusion (8.5) with a = 0 and constant diffusivity 
c(x) in a two-dimensional channel (8.4). To be specific, we set L = 5, // = 1, and 
the random bottom boundary as a random field with zero mean and an exponential 
two-point covariance function 

0 < r, s < L, (8.8) 

where b > 0 is the correlation length. In the computational examples below, b = 
L/5 = 1, which corresponds to a boundary of moderate roughness. Finally, we 
prescribe Dirichlet boundary conditions u = I at X 2 = 11 and u — 0 elsewhere. 


Chh{r, s) = E [h(r, co)h(s, co)] = a 2 exp 
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(a) (b) 

Figure 8.1 Channels with a rough wall generated with the 10-term ( K = 10) KL expan- 
sion (8.9). (a) Four sample realizations of the bottom boundary s(xi, a>j) (j = 
1 , .... 4). (b) A sample realization of the channel in the physical domain (xi , x 2 ) 
and in the mapped domain in (£ l5 £ 2 ). Chebyshev meshes are used in both do- 
mains. (More details are in [130].) 


We employ the finite-term Karhunen-Loeve (KL)-type expansion (4.8) to 
decompose the boundary process. That is, 

d 

h(x\, co) a E y/h'l'kixi )Z k (ca), (8.9) 

k=i 

where {A,*, t//>(xi)} are the eigenvalues and eigenfunctions of the integral equations 

L 

Chh(r, s)-^k(r)dr — A^i/rtC?), k=\ (8.10) 

We further set { Z, (o>)} ~ (/(— 1 , 1 ) to be independent uniform random variables in 
(—1, 1) and use the parameter 0 < a < 1 to control the maximum deviation of the 
randomness. (In the computational examples in this section, we set a = 0.1.) We 
employ Legendre polynomials as the gPC basis functions. 

It is worthwhile to stress again that the expansion (8.9) introduces two sources of 
errors — errors due to the finite d - term truncation and errors due to the assumption 
of independence of {Za-(<w)}. The truncation error is typically controlled by select- 
ing the value of d to ensure that the eigenvalues {A^.} with k > d are sufficiently 
small. For example, in this example the expansion with d = 10 captures 95 percent 
of the total spectrum (i.e., eigenvalues). A few realizations of the bottom boundary, 
obtained by the 10-term KL expansion, are shown in figure 8.1a. In figure 8.1b, 
one realization of the channel domain is mapped onto the corresponding rectangu- 
lar domain E — [0, L] x [0, H], Also shown here are the Chebyshev collocation 
mesh points that are used to solve the mapped stochastic diffusion problem (8.6). 

Figure 8.2 shows the first two moments of the solution, i.e., its mean (top) and 
standard deviation (STD) (bottom). We observe that the STD reaches its maximum 
close to, but not at, the random bottom boundary. 
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Figure 8.2 The mean and the STD of the dependent variable u computed with the stochastic 
Galerkin method. 


To ascertain the convergence of the polynomial chaos expansion, we examine the 
STD profile along the cross section y = 0.25, where the STD is close to its maxi- 
mum. Figure 8.3 shows the STD profiles obtained with different orders of Legendre 
expansions. One can see that the second order is sufficient for the Legendre expan- 
sion to converge. Although not shown here, the convergence of the mean solution 
is similar to that of the STD. 

Monte Carlo simulations (MCSs) are also conducted to verify the results ob- 
tained by the stochastic Galerkin method. Figure 8.4 compares the STD profile 
along the cross section y — 0.25 computed via the second-order Legendre expan- 
sion with those obtained from MCS. We observe that as the number of realizations 
increases, the MCS results converge to the converged SG results. With about 2,000 
realizations, the MCS results agree well with the SG results. In this case, at second 
order ( N — 2) and 10 random dimensions id = 10), the gPC stochastic Galerkin 
method requires ( N + d)\/N\d\ = 66 basis functions and is computationally more 
efficient than Monte Carlo simulation. 

Often analytical mapping is not available; then a numerical technique can be em- 
ployed to determine the mapping ([130]). This involves solving of a set of boundary 
value problems for the mapping. Other techniques for casting random domain prob- 
lems into deterministic problems include the boundary perturbation method [126], 
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Figure 8.3 The STD profiles along the cross section y = 0.25 computed with the first-, 
second-, and third-degree Legendre polynomials. 



Figure 8.4 The STD profiles along the cross section y = 0.25 computed with the SG method 
(second-order Legendre chaos) and Monte MCS consisting of 100, 500, 1000, 
and 2000 realizations. 
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isoparametric mapping [15], the fictitious domain method [14], the extended finite 
element method [82], and a Lagrangian approach that works well for solid defor- 
mation [2]. Interested readers are strongly encouraged to consult the references. 


8.2 BAYESIAN INVERSE APPROACH FOR PARAMETER ESTIMATION 

When solving a stochastic system, e.g., (4.13), it is important that the probability 
distribution functions of the input random variables Z be available. This, along 
with the independence condition, allows us to sample the inputs or to build basis 
functions (such as gPC) to solve the system. A question naturally arises: What if 
there is not enough available information to determine and specify the parameter 
distributions? This is a problem of practical concern because in many cases there 
are not enough measurement data on the parameters — some parameters cannot even 
be measured. However, sometimes other measurement data are available — data not 
on the parameters but on some other quantities that can be computed. In this case, 
an inverse parameter estimation can be carried out to estimate the true distributions 
of the input random parameters. 

The field of parameter estimation is not new — much research has gone into it 
for decades. Here, however, we will present an approach using Bayesian inference. 
More importantly, the approach is built upon the gPC algorithms in such way that 
it does not incur any simulation effort in addition to a one-time gPC simulation. 
In other words, for the inverse problem here, a forward gPC simulation is required 
only once, and the rest becomes off-line postprocessing. The strong approximation 
property of gPC expansion plays a vital role here. 

Let us assume that each random variable Z, has a prior distribution F, (z, ) = 
P(Zi < Zi) e [0, 1]. The distribution can be made upon assumption, intuition, or 
even speculation when there are not sufficient data. Here we focus on continu- 
ous random variables. Subsequently, each Z, has a probability density function 
jTj(zi) = d F t ( z, ) /d z, . We also assume the variables are mutually independent — 
another assumption when not enough data are available to suggest otherwise. Thus, 
the joint prior density function for Z is 

n z {z) = Y\ni{zi). (8.11) 

i=t 

Whenever possible, we will neglect the subscript of each probability density and 
use 7r(z) to denote the probability density function of the random variable Z, rc/(z), 
unless confusion would arise. 

Let 

d‘ = g(u)eW d (8.12) 

be a set of variables that one observes, where g : R"“ —> W' d is a function re- 
lating the solution u to the true observable d' . We then define a forward model 
G : R"* — »■ R"“' to describe the relation between the random parameters Z and the 
observable d ‘ : 


d' = G(Z ) = go u(Z). 


(8.13) 
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In practice, measurement error is inevitable and the observed data d may not match 
the true value of d' . Assuming additive observational errors, we have 

d = d‘+e= G(Z) + e, (8.14) 


where e e K' d are mutually independent random variables with probability den- 
sity functions n(e) = rT/=i 71 ( e i)- We make the usual assumption that e are also 
independent of Z. 

The Bayesian approach seeks to estimate the parameters Z when given a set of 
observations d. To this end, Bayes’ rule takes the form 


jt(z\d) = 


jt{d\z)n{z) 
f Tt{d\z)Tt(z)dz 


(8.15) 


where j r(z) is the prior probability density of Z; n{d\z) is the likelihood function; 
and n{z\d), the density of Z conditioned on the data d, is the posterior probability 
density of Z. For notational convenience, we will use n d (z) to denote the posterior 
density jt(z\d) and L{z) to denote the likelihood function ic(d\z). That is, (8.15) 
can be written as 


7T d (z) 


L(z)n(z) 

J L(z)7t(z)dz 


(8.16) 


Following the independence assumption on the measurement noise e, the likelihood 
function is 


nd 

L{z) = Ji(d\z) = ]~~[ n ei id, - Gj(z)). (8.17) 


The formulation, albeit simple, poses a challenge in practice. This is largely be- 
cause the posterior distribution Tt d does not have a closed and explicit form, thus 
preventing one from sampling it directly. A large amount of literature has been 
devoted to this challenge, with one of the most widely used approaches being the 
Markov chain Monte Carlo method (MCMC). For an extensive review, see [101]. 
Since most of the approaches are based on sampling, the main concern is to im- 
prove efficiency because each sampling point requires a solution of the underlying 
forward problem G(Z) and can be time-consuming. 

The gPC-based numerical method, in addition to its efficiency in solving the 
forward problem, provides another remarkable advantage here. For all 1 < i < nj, 
let G Nj(Z) be a gPC approximation for the / th component of G(Z), 


N 

GnAZ)= ^a k>i O k (Z), 

|k|=0 


(8.18) 


where the expansion coefficients {a k ,, } are obtained by either a stochastic Galerkin 
or a stochastic collocation method. Then we effectively have an analytical repre- 
sentation of the forward problem in terms of Z, which can be sampled at an arbi- 
trarily large number of nodes by simply evaluating the polynomial expression at the 
nodes. Thus, if we substitute the gPC approximation into the likelihood function, 
we obtain an approximate posterior density that can be easily sampled without any 
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simulation effort, in addition to the stochastic forward problem which needs to be 
solved only once. 

The gPC approximation of the posterior probability is 


n d N (z) = 


L n (z)jt(z ) 

/ LN(z)7T{z)dz’ 


(8.19) 


where it(z) is again the prior density of Z and is the approximate likelihood 
function defined as 


nj 

L n (z) 4 n N {d\z) = Y\ji ei (di - G nj (z)). (8.20) 


The error of the approximation can be quantified by using Kullback-Leibler di- 
vergence (KLD), which measures the difference between probability distributions 
and is defined, for probability density functions jt\ (z) and 7r 2 (z), as 

D(txi\\tx 2 )= f Jti(z) log T ^-y^-dz. (8.21) 

J 7T 2 (z) 

It is always nonnegative, and D(jt\ ||7r 2 ) = 0 when jt\ = 7T 2 . 

Under the common assumption that the observational error in (8.14) is indepen- 
dently and identically distributed (i.i.d.) Gaussian, e.g., 

e ~ N(0, cr 2 I), (8.22) 


where o > 0 is the standard deviation and I is an identity matrix of size «,/ x nj, 
the following results can be established. 


Lemma 8.4. Assume that the observational error in (8.14) has an i.i.d. Gaussian 
distribution (8.22). If the gPC expansion Gy in (8.18) of the forward model con- 
verges to Gi, 

\\G i (Z)-G N i(Z)\\ L 2 -*0, 1 <i<n d , N -» 00 , 

’ n Z 

then the posterior probability Tt d N in (8.19) converges to the true posterior probabil- 
ity jt d in (8.16) in the sense that the Kullback-Leibler divergence (8.21) converges 
to zero; i.e., 

D(n d N \\Tt d ) -* 0, N 00 . (8.23) 

Theorem 8.5. Assume that the observational error in (8.14) has an i.i.d. Gaussian 
distribution (8.22) and that the gPC expansion G^j in (8.18) of the forward model 
converges to Gi in the form of 

\\Gi(Z)-G Nj (Z)\\ L 2 < CN~ a , 1 <i<n d , 

TC Z 

where C is a constan t independent of N and a > 0, then for sufficiently large N, 
D(n d N \\7t d ) < N~ a . (8.24) 

Proofs of these results can be found in [76]. 
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Figure 8.5 The posterior distribution of the boundary condition uncertainty <5 (denoted as 
p(8\data) here) of the Burgers’ equation. The exact posterior from Bayes’ rule 
is shown here along with the numerical results from the gPC collocation-based 
algorithms with order N = 4 and N = 8. Also shown is the true (and unknown) 
perturbation of 8'. 


Example 8.6 (Burgers’ equation). Let us return to the Burgers’ equation exam- 
ple in section 1.1.1. As shown in figure 1.1, a small amount of uncertainty with 
S ~ U((), 0.1) can produce a large response in the location of the transition layer 
because of the supersensitive nature of the problem. On the other hand, the pre- 
diction produces the distribution range of the output, albeit correctly, too big to be 
useful. This suggests that the assumption about the uncertainty of S is too wide and 
should be refined both in its range and distribution. This can be achieved using the 
Bayesian inverse estimation when some observation data are available. 

Without actual experimental data, we generate “data” numerically. This is ac- 
complished by fixing a “true” perturbation S‘ , conducting a high-order determinis- 
tic simulation to compute the true location of the transition layer d 1 at steady state, 
and then adding Gaussian noise e to produce the data d = d‘ + e. The data d are 
then used for the Bayesian inverse estimate for the posterior distribution of 8. 

In figure 8.5, the numerical results of the posterior density are shown with gPC 
orders of N = 4 and N = 8. Compared to the exact posterior density obtained 
from the Bayesian rule directly (in this case the exact solution can be solved), we 
notice that the gPC results converge. For reference, the true and yet unknown lo- 
cation of the perturbation S' is also plotted. We observe that the posterior density 
clusters around the truth, as expected. The convergence of the gPC -based Bayesian 
algorithm is examined in more detail in figure 8.6 for orders as high as N — 200. 
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N 


Figure 8.6 Convergence of the numerical posterior distribution by gPC-based Bayesian al- 
gorithms at increasing orders, along with the convergence of the gPC forward 
model prediction. 


We observe the familiar exponential convergence. The convergence of the gPC for- 
ward problem solver is also plotted. It is clear that the posterior density converges 
at least as fast as (in fact, even faster than) the forward problem. This is consistent 
with the theoretical result. For more details on the analysis and numerical results, 
see [76]. 


8.3 DATA ASSIMILATION BY THE ENSEMBLE KALMAN FILTER 

Any mathematical and numerical models, deterministic or stochastic, no matter 
how sophisticated, are approximations to the true physics. Though many models 
are accurate for a wide range of spatial and time domains, many can deviate from 
the “truth” quickly. (For example, weather forecasting.) In addition to improving 
the models, another way to improve the prediction is to take advantage of measure- 
ment data, which reflect the physical truth, sometimes partially (scarce and indirect 
measurement) and approximately (inaccurate measurement). 

In data assimilation, data arrive sequentially in time, and one seeks to incor- 
porate both the data and the prediction of the mathematical/numerical models to 
produce better predictions. There has been an extremely large amount of litera- 
ture on various methods of data assimilation, with the most popular ones based on 
either a Kalman filter (KF) [54] or a particle filter. Here we discuss only the ensem- 
ble Kalman filter (EnKF) [28], a variance of the Kalman filter, and explain how 
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gPC-based stochastic methods can be used to significantly improve the perfor- 
mance of the EnKF. The notation here is somewhat different from that in the rest 
of the book, particularly the use of bold letters, as we try to follow the notations 
commonly used in the data assimilation literature. 

Let u ^ e R"' , m > 1 , be a vector of forecast state variables (denoted by the 
superscript /) that are modeled by the following system: 


— = /(fV), fe(OJ], (8.25) 

at 

U R(O) = u 0 (Z), (8.26) 


where T > 0 and Z e S, d , d > 1, is a set of random variables parameterizing the 
random initial condition. The model (8.25) and (8.26) is obviously not a perfect 
model for the true physics, and the forecast may not represent the true state vari- 
ables, u' e R m , sufficiently well. If a set of measurements d e R*, i > 1, are 
available as 


d = Hu' + e, (8.27) 

where H : W — » R f is a measurement operator relating the true state variables 
u' and the observation vector d e R' and e e R f> is measurement error. Note 
that the measurement operator can be nonlinear, although it is written here in a 
linear fashion by following the traditional exposition of the (ensemble) Kalman 
filter. Also, characterization of the true state variables u' can be highly nontrivial in 
practice. Here we assume that they are well-defined variables with dimension m . 

The objective of data assimilation is to construct an optimal estimate of the true 
state, the analyzed state vector denoted as u° e R" ! , based on the forecast u 1 and 
the observation d. Note that it is possible to add a noise term in (8.25) as a model 
for the modeling error. Here we restrict ourselves to the deterministic model (8.25) 
with the random initial condition (8.26). 


8.3.1 The Kalman Filter and the Ensemble Kalman Filter 

The Kalman filter is a sequential data assimilation method that consists of two 
stages at each time level when data are available — a forecast stage where the system 
(8.25) and (8.26) is solved and an analysis stage where the analyzed state u" is 
obtained. Let pR e R" iX "' be the covariance matrix of the forecast solution u ' . 
The analyzed solution u" in the standard KF is determined as a combination of the 
forecast solution uR and the measurement d in the following manner: 

u a =uR + K(d-H U R), (8.28) 

where K is the Kalman gain matrix defined as 

K = pRH r (HP / H r + Rr 1 . (8.29) 

Here the superscript T denotes matrix transpose, and R e R ,x< is the covariance 
of the measurement error e . The covariance function of the analyzed state u". P" e 
M mxm , is then obtained by 

P a = (I - KH)PR (I - KH) r + KRK r = (I - KH)pR, 


(8.30) 
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where I is the identity matrix. When the system (8.25) is linear, the KF can be 
applied in a straightforward manner, as equations for evolution of the solution co- 
variance can be derived. For nonlinear systems, explicit derivation of the equations 
for the covariance function is not possible. Subsequently, various approximations 
such as the extended Kalman filter (EKF) were developed. Their applicability is 
limited to a certain degree depending on the approximation procedure. Further- 
more, in practical applications, forwarding the covariance functions (8.30) in time 
requires explicit storage and computation of P 7 , which scales as 0(m 2 ) and can be 
inefficient when the dimension of the model states, m, is large. 

The ensemble Kalman filter (EnKF) ([11, 28]) overcomes the limitations of the 
Kalman filter by using an ensemble approximation of the random state solutions. 
Let 

(u 7 );, i = l,...,M, M> 1, (8.31) 

be an ensemble of the forecast state variables u 7 , where each ensemble member is 
indexed by the subscript i = I .... . M and is obtained by solving the full nonlinear 
system (8.25). The analysis step for the EnKF consists of the following update 
performed on each of the model state ensemble members 

(u a \ = (u 7 ), + K e ((d),- - H(u 7 ),), i = 1 M, (8.32) 


where 

K e = P 7 H r (HP 7 H r + Rg)-' (8.33) 


is the ensemble Kalman gain matrix. Here 

p 7 4 (u7 - n 7 )(u/ - u 7 y ~ p 7 , 

P° e = (U a -U°)(U° -VL a y ~ P fl , 


(8.34) 


are the approximate forecast covariance and analysis covariance, respectively, 
obtained by using statistical averages of the solution ensemble (denoted by the 
overbar), and R e = ee T — R is the approximate observation error covariance. 
Therefore, the covariance functions are approximated by ensemble averages and 
do not need to be forwarded in time explicitly. An extensive review of the EnKF 
can be found in [29]. 


8.3.2 Error Bound of the EnKF 

As an approximation to the Kalman filter, one obvious source of error for the EnKF 
is from the sampling. Note that here we define error as numerical error — it is the 
difference between the result obtained by the KF and that obtained by the EnKF. 
This is different from modeling error between the result of the KF and the physi- 
cal truth. For EnKF, the overall error consists of both the numerical error and the 
modeling error. Here we study only the former. 

To understand the impact of numerical error more precisely, we cite here an error 
bound of the EnKF derived in [66]. Let t\ < ti < ■ ■ ■ be discrete time instances 
at which data arrive sequentially and assimilation is made. Without loss of gen- 
erality, let us assume that they are uniformly distributed with a constant step size 
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AT = tk — tk~i, Vk > 1. Let E„ be the numerical error of the EnKF, that is, the 
difference between the EnKF results and the exact KF results measured in a proper 
norm at time level t n , n > 1, then the following bound holds: 

E n < ^ e^j exp (A -t„), (8.35) 

where Eq is the error of sampling the initial state, e k is the local error at time level 
tk, 1 < k < n, and A > 0 is a constant. The local error scales as 

e k ~ O (At p , aM ~ a ) , At -* 0, M -» 00 , (8.36) 

where 0(At p ) denotes the numerical integration error in time induced by solving 
(8.25) and (8.26) with a time step At and a temporal integration order p > 1, 
a > 0 is the noise level of the measurement data and scales with the standard 
deviation of the measurement noise, M is the size of the ensemble, and a > 0 is 
the convergence rate of the sampling scheme. For Monte Carlo sampling, a — 1/2. 
In most cases, this sampling error dominates. A notable result is that the constant A 
depends on the size of the assimilation step in an inverse manner, i.e., A ot AT~ l . 
This implies that more frequent data assimilation by the EnKF can magnify the 
numerical errors. Since more frequent assimilation is always desirable (whenever 
data are available) for a better estimate of the true state, it is imperative to keep 
the numerical errors, particularly the sampling errors, of the EnKF under control. 
Although the sampling errors can be easily reduced by increasing the ensemble 
size, in practice this can significantly increase the computational burden, especially 
for large-scale problems. 


8.3.3 Improved EnKF via gPC Methods 


Here we demonstrate that one can yet again take advantage of a highly accurate 
gPC approximation to construct an improved EnKF scheme with much reduced 
numerical error. Let 

N 

u ^(t,Z) = ^fi i / (f)«l» i (Z) (8.37) 

|i|=0 


be the gPC solution to the forecast equations (8.25) and (8.26) with sufficiently 
high accuracy, where the expansion coefficients u/ (f) can be obtained by either the 
gPC Galerkin procedure or the gPC collocation procedure. 

In addition to offering efficiency for the forecast solution, another (often over- 
looked) advantage of gPC expansion is that it provides an analytical representation 
of the solution in terms of the random inputs. All statistical information about u(, 
can be obtained analytically or with minimum computational effort. For example, 
the mean and covariance are 


57/ 


= £ 



(8.38) 


respectively, where y j = E[0?]. And they can be used as accurate approximations 
of the exact mean and covariance of the forecast solution u L Furthermore, one can 
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generate an ensemble of solution realizations by sampling the random variables 
Z in (8.37). This procedure involves nothing but polynomial evaluations, and thus 
generating an ensemble with an arbitrarily large number of samples does not require 
any computation of the original governing equations (8.25) and (8.26). Let 

N 

(u£),- = Y, Uk(04» k ((2)»), i = M » 1, (8.39) 

|k|=0 

be an ensemble of the forecast solution realizations with size M, where (Z), , i = 
1, . . . , M, are Monte Carlo samples of the random vector Z. Equipped with a 
knowledge of the solution statistics, particularly the mean and covariance from 
(8.38), we can apply the EnKF scheme (8.32) to obtain analyzed states. 

(u fl w ); = (Uy), + K ;V (d, - H(u£)i), i = \, ... ,M, (8.40) 

where K v is the gPC Kalman gain matrix defined as 

K.V = P^H r (HP^H 7 " + R) _1 , (8.41) 

which approximates the Kalman gain matrix (8.29). 

The key ingredients and advantages of the gPC-based EnKF are 

• Solution of the forecast problem by a gPC-based method, either a Galerkin or 
a collocation method. This, in many cases, is (much) more efficient than the 
traditional Monte Carlo sampling employed in the EnKF. 

• At the update stage, one can use the gPC forward problem solution to generate 
an arbitrarily large number of samples and update them individually by (8.40), 
similar to the traditional EnKF. The ability to update a large ensemble of solu- 
tions results in a significant reduction in sampling errors. Moreover, this step 
is virtually “free,” in the sense that generating the large solution ensemble via 
gPC is nothing but evaluation of the gPC polynomial expression and induces 
no simulation cost. 

• For a linear system of equations (8.25) with Gaussian noise, the gPC-based 
EnKF becomes equivalent to the standard Kalman filter. 

More details about these methods can be found in [67], 
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Appendix A 


Some Important Orthogonal Polynomials 
in the Askey Scheme 


Here we summarize the definitions and properties of some important orthogonal 
polynomials from the Askey scheme. Denote { Q n (x)} as an orthogonal polynomial 
system with the orthogonal relation 

J Q n (x)Q m (x)w(x)dx = h„8 m „ 

for continuous x, or in the discrete case, 


^2 Qn(x)Qm(x)w(x) = 


where S is the support of w(x). The three-term recurrence relation takes the form 
x Qn (x) = b„Q n+ i(x) + y„Q„(x) + c„Q„- i(x), n > 0, 
with initial conditions i(x) = 0 and Qq(x) = 1. Another way of expressing the 
recurrence relation is 


Q„+i(x) = (A„x + B„)Q„(x) - C„Q„-i(x), n > 0, (A.l) 

where A„, C„ f 0 and C n A n A„_\ > 0. It is straightforward to show that, if we 
scale variable x by denoting y — ax for a > 0, then the recurrence relation takes 
the form 

Sm-iOO = {A n y + aB n )S n {y)-ot 2 C n S n _Xy)- (A.2) 


Another important property is that these orthogonal polynomials are solutions of a 
differential equation 


s(x)y" + r(x)y + Xy — 0 (A. 3) 

in continuous cases, and a difference equation 

i(x)AVy(x) + r(x)A_y(x) + Xy(x) = 0 (A.4) 

in discrete cases, where s(x) and r(x) are polynomials of at most second and first 
degree, respectively, and X is a constant. The notations for the discrete cases are 

A fix) = f{x + 1) - fix), V fix) = fix) - fix - 1). 


When 


X = X n = — nr ' 


1 

-nin — 1)5 , 
2 


the equations have a particular solution of the form _y(x) = Q„ix), which is a 
polynomial of degree n . 
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A.l CONTINUOUS POLYNOMIALS 


A. 1.1 Hermite Polynomial //„ (x) and Gaussian Distribution 

Definition: 


H„(x) = (2x) n 2 F 0 




Orthogonality: 


where 


Recurrence relation: 



H m (x) H n (x)w(x)dx = n\S m „, 


1 

w(x ) = e 

■Jin 


x 2 /2 


H n+ i(x) = xH n (x) — nH„-\ (x). 

Rodriguez formula: 

e- x2 ' 2 H n (x) = (-1)"^ (e“* 2/2 ) . 

Differential equation: 

y”(x) — x)/ (x) + ny(x) = 0, y(x) = H n (x). 


(A. 5) 

(A. 6) 

(A. 7) 

(A. 8) 

(A. 9) 

(A. 10) 


A. 1.2 Laguerre Polynomial L^ a) (x) and Gamma Distribution 

Definition: 

L' a> (x) = - — — — — i Fi(—n; a + 1; x). 
n\ 


Orthogonality: 

L ^ (x ) L { “Hx)w(x)dx 

where 



(«+!)» 


a > — 1, 


w(x) = 


x a e x 

r(a+ 1)' 


Recurrence relation: 

(« + l)L^j(x) — (2 n + a + 1 — x)L ( “\x) + (n + a)L\®\(x) = 0. 
Normalized recurrence relation: 

xq n {x) = q n+ \(x) + (2 n + a + \)q„(x) + n(n + u)q n - i(x), 


where 


Z,< o) (jc) 


(- 1 )' 


q n (x). 


(A.l 1) 

(A. 12) 

(A. 13) 

(A. 14) 

(A. 15) 


n ! 
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Rodriguez formula: 


e X x a L [ f\x) = (e A x” + “) 

n\ ax n 


(A. 16) 


Differential equation: 

xy"{x) + (a + 1 — x)/(x) + ny(x) = 0, y(x) = L^fx). (A. 17) 
Recall that the gamma distribution has the probability density function 


fix) = 


„-x/P 


a > — 1, B > 0. 


(A. 18) 


/!“+' T (a + 1)’ 

The weighting function of Laguerre polynomial (A. 13) is the same as that of the 
gamma distribution with the scale parameter (3 = 1 . 


A.1.3 Jacobi Polynomial Pf’^fx) and Beta Distribution 

Definition: 

(x) = 

n\ 

Orthogonality: 


F\ I —n ,n a + ft + 1 ; a + 1; 


1 — x 


J P^’ w (x) J P„ ( “’ ffl (x)tn(x)c/x = h 2 n S mn , a, P > -1, 


where 


hi = 


w(x) = 


Recurrence relation: 


(a + 1)„ (£ + !)„ 


n\(2n + a + P + l)(a + P + 2)„_i ’ 

r(a + p + 2) 


2“+^+ 1 r(a + l)T(yS + 1) 


xP^ix) = — + 


(2« + a + p + 1 ) (2/7 + a + p + 2) 
P 2 -a 2 


(2 n + a + P){2n + a + P + 2) 
2(n + a)(n + P) 

(2 n + a + P)(2n + a + P + 1) 
Normalized recurrence relation: 

«2 _ 


P^\x) 


P^\x). 


XPn (x) = p n+ i(x) + 


— a 


(2 n + a + P)(2n + a + P + 2) 
4n(« + a)(« + P)(n + a + P) 


Pn(x) 


(A. 19) 


(A. 20) 


(1 -x)“(l + x)f (A.21) 


(A. 22) 


(2 n + a + ft — \)(2n + a + ft) 2 (2n + a + ft + 1) 


ft-tW, (A. 23) 


r , ( , w = 2L±L+ALiL RW . 

2"n! 


where 
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Rodriguez formula: 

(1 - JC)“(1 + xfp^ix) = [d - x)” + “(l + x) n+fi ] . (A.24) 

2 n n\ dx n 

Differential equation: 

(1— x 2 )y"(x) + [/0 — a — (a + P+2)x]y'(x)+n(n+oi + P+l)y(x) — 0, (A.25) 

where _y(x) = P^ a '^\x). 


A.2 DISCRETE POLYNOMIALS 


A.2.1 Charlier Polynomial C„(x; a) and Poisson Distribution 

Definition: 

/ 1 

C n (x; a) — 2 F 0 I — n, —x; ; 

V a 

Orthogonality: 


(A. 26) 


'S~' — C m (x; a)C„{x\ a) — a "e a n\S m „, a > 0. (A. 27) 

' x\ 

x=0 

Recurrence relation: 

—xC n (x; a) = aC„ + i(x; a) — (n + a)C„(x; a) + wC„_i(x; a). (A. 28) 

Rodriguez formula: 

a x „ ( a x \ 

— -C„(x; a) = V ( — r ) . (A. 29) 

x! \x! / 

Difference equation: 

— ny(x ) = ay(x + 1) — (x + a)y(x ) + xy(x — 1), y(x) — C„(x; a). (A. 30) 

The probability function of Poisson distribution is 

f(x\a) = e~ a — , k = 0, 1,2, .... (A.31) 

x! 

Despite a constant factor e~ a , it is the same as the weighting function of Charlier 
polynomials. 


A.2.2 Krawtchouk Polynomial K lt (x; p, N ) and Binomial Distribution 

Definition: 


K n (x\ p, N ) 


2 Li 



n = 0, (A. 32) 
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Orthogonality: 


N , , 

Y ( Y )^ X(1 - P) N ~ x K m {x\ P, N)K n (x ; p, N ) 


(— 1 )”«! f 1 — p 
(~N)„ 


S mn , 0 < p < 1. 


(A. 33) 


Recurrence relation: 


—xK(x; p , N) = p(N -n)K„+i(x- p , N) — [p(N — n) + n(\ — p)]K n {x\ p , N) 
+ »(1 — p)K n ^\{x\ p, N). (A. 34) 


Rodriguez formula: 




K n (x\ p, N ) = V" 


Difference equation: 



(A.35) 


— ny{x) — p(N — x)y(x + 1) — [p(N — x) — xq\y(x) + xqy(x — 1), (A. 36) 


where y(x) = K„(x; p, N ) and q = 1 — p. 

Clearly, the weighting function from (A. 33) is the probability function of the 
binomial distribution. 


A.2.3 Meixner Polynomial M„( x; (i, c ) and Negative Binomial Distribution 

Definition: 


M n {x\ p, c) 



(A.37) 


Orthogonality: 


(P)x v C~ n n 

Y — —c M m (x\ P, c)M n (x; p, c) = 


03)h(1 — c)P 


>0, 0 < c < 1. 

(A.38) 


Recurrence relation: 


(c — 1 )xM n (x; P, c ) = c(n + P)M n+ \{x\ ft, c ) — [n + (n + P)c]M n (x\ ft, c) 

+ nM„_ l (x; /3,c). (A.39) 


Rodriguez formula: 


(P)xC X 

x\ 


M n (x\ P,c) = V" 


(P + n) x c x 
x! 


(A. 40) 


Difference equation: 

n{c - 1 )y(x) = c(x + P)y(x + 1) - [x + (x + P)c]y(x) + xy(x - 1), (A.41) 
where y{x) — M n (x\ P, c). 
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The weighting function is 

/(*)= ^(l-c)V, 0 < c < 1 , /? > 0 , x = 0 , 1 , 2 ,.... (A. 42) 

x! 

It can be verified that it is the probability function of a negative binomial distribu- 
tion. In the case where /J is an integer, it is often called the Pascal distribution. 

A.2.4 Hahn Polynomial Q n (x; a, fi, N) and Hypergeometric Distribution 

Definition: 

Q n (x; a , /3, A) = 3 T))— «, n + a + /i + 1, — x; a+ 1, — A ; 1), n = 0, 1, . . . , A. 

(A. 43) 

Orthogonality: Fora > — 1 and fi > — 1 or for a < — A and fi < —A, 

f a + x\ f 6 + N — x\ 0 

x J (j v _ x JQ m (x-a,p,N)Q„(x-a^,N) = h 2 n S mn , (A.44) 

where 

^2 _ ( — 1 )" ( n + a + P + 1)jv+i(^ + l)« w ! 

7,1 ~~ (2n +a + /? + l)(a + 1)„(-A)„A! ' 

Recurrence relation: 

-xft(x) = A n Q n+ i(x) - (4„ + C n )Q„(x) + C,g,. 1 (x), (A.45) 

where 


6»(x) := Q n (x; a, fi, A) 


and 

A„ 

C„ 

Rodriguez formula: 


(« + a + + 1 )(n + a + 1)(A — n) 

(2 n + a + fi + l)(2w + a + fi + 2) 

«(« + a + fi + A + l)(w + fi) 

(2n + a + fi)(2n + a + fi + 1) 


w(x; a, /8, A) g„(x; a, /J, A) 
where 


(-1)"03 + 1), 
(-#)» 


V"[hj(x; a + n, fi + n, N 


- «)L 
(A. 46) 


ut(x; a, /S, A) 


fP + A - x\ 

V / 


Difference equation: 


n(n + a + fi + 1 )y(x) = B(x)y(x + 1) — [ B(x ) + D(x)]y(x) + D(x)y(x — 1), 

(A. 47) 


where y{x) = Q„(x; a, fi , A), B(x ) = (x + a + l)(x — A), and D(x) = x(x — 

P — A*— 1). 
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If we set a = —a — 1 and — —fi — 1, we obtain 


ui(x) = 


!)(A) 


(N-a-0- 1\ 


1+P\ 


Apart from the constant factor I /( v “ v ^ this is the definition of hyper- 
geometric distribution. 
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Appendix B 


The Truncated Gaussian Model G(a, /3) 


The truncated Gaussian model was developed in [123] in order to circumvent the 
mathematical difficulty resulting from the tails of Gaussian distribution. It is an 
approximation of Gaussian distribution by generalized polynomial chaos (gPC) Ja- 
cobi expansion. The approximation can be improved either by increasing the order 
of the gPC expansion or by adjusting the parameters in the Jacobi polynomials. 
The important property of the model is that it has bounded support, i.e., no tails. 
This can be used as an alternative in practical applications where the random inputs 
resemble Gaussian distribution and the boundedness of the supports is critical to 
the solution procedure. 

While the procedure of approximating (weakly) a Gaussian distribution by Ja- 
cobi polynomials was explained in section 5.1.2. here we tabulate the results for fu- 
ture reference. The gPC Jacobi approximation for Gaussian distribution is denoted 
as G(a, P) with a, ft > —1. Because of the symmetry of Gaussian distribution, we 
set a — p in the Jacobi polynomials. 

In figures B.1-B.3, the probability density functions (PDFs) of the gPC Jacobi 
chaos approximations are plotted for values of a = ft = 0 to 10. For a = (3 = 0, 
Jacobi chaos becomes Legendre chaos, and the first-order expansion is simply a 
uniform random variable. In this case, Gibbs’ oscillations are observed. As the 
values of (a, ft) increase, the approximations improve. The expansion coefficients 
at different orders are tabulated in table B.l, together with the errors in variance 
and kurtosis compared with the exact Gaussian distribution. It is seen that, with 



Figure B.l Gaussian random variables approximated by Jacobi chaos. Left: a = ft = 0. 
Right: a = ft = 2. 
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Figure B.2 Gaussian random variables approximated by Jacobi chaos. Left: a = f) — 4. 
Right: a = /S = 6. 



Figure B.3 Gaussian random variables approximated by Jacobi chaos. Left: a = f) — 8. 
Right: a = fi = 10. 


a = ji = 10, even the first-order approximation, which is simply a beta random 
variable, has an error in variance of as little as 0. 1 percent. The errors in kurtosis are 
larger because the Jacobi chaos approximations do not possess tails. This, however, 
is exactly our objective. 


Table B.l Approximating Gaussian Random Variables via Jacobi Chaos: Expansion Coefficients yk and errors" 



a = f3 = 0 

C4 

II 

SA 

II 

a 

■St- 

II 

'Q. 

II 

a 

a = (3 — 6 

a = P = 8 

0 

II 

<3Q. 

II 

a 

Vi 

1.69248 

8.78271-1) 

6.62181-1) 

5.52731-1) 

4.83991-1) 

4.35751-1) 

C2 

4.51704(— 2) 

8.253461—3) 

3.463011—3) 

2.007291-3) 

1.388421-3) 

1.072311-3) 

C4 

1.35894 

7.050241-1) 

4.790891-1) 

3.635571—1) 

2.932461-1) 

2.459161-1) 

.V3 

4.83991— 1) 

7.54931-2) 

2.60111-2) 

1.22161-2) 

6.779701-3) 

4.177921-3) 

C2 

1.17071(-2) 

8.518161-4) 

4.492451-4) 

4.239831-4) 

4.338941-4) 

4.452821-4) 

C4 

5.020971-1) 

7.97474(-2) 

3.33201(-2) 

2.400641-2) 

2.21484(-2) 

2.225391-2) 

V5 

2.70641-1) 

1.99591-2) 

2.99361-3) 

2.3531(-4) 

—3.308881—4) 

-4.195391-4) 

C2 

5.048381-3) 

3.970591-4) 

3.968801-4) 

4.229031-4) 

4.282831-4) 

4.250431-4) 

<?4 

2.555261-1) 

2.293731-2) 

1.92101(-2) 

2.150951-2) 

2.068461-2) 

2.083171-2) 


“62 is the error in variance; e 4 is the error in kurtosis. There is no error in the mean, yt = 0 when k is even. 
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supersensitivity. See Burgers’ equation 

three-term recurrence, 26, 105 

Vandermonde matrix, 37-38, 80 

Weierstrass theorem, 30 


